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ABSTRACT 



Aims. We investigate the building of unified models that can predict the matter-density power spectrum and the two- 
point correlation function from very large to small scales, being consistent with perturbation theory at low k and with 
halo models at high k. 

Methods. We use a Lagrangian framework to re-interpret the halo model and to decompose the power spectrum into 
"2-halo" and "1-halo" contributions, related to "perturbative" and "non-perturbative" terms. We describe a simple 
implementation of this model and present a detailed comparison with numerical simulations, from k ~ 0.02 up to 
100/iMpc" 1 , and from x ~ 0.02 up to 150/i _1 Mpc. 

Results. We show that the 1-halo contribution contains a counterterm that ensures a k 2 tail at low k and is important 
not to spoil the predictions on the scales probed by baryon acoustic oscillations, k ~ 0.02 to 0.3/iMpc - . On the 
other hand, we show that standard perturbation theory is inadequate for the 2-halo contribution, because higher order 
terms grow too fast at high k, so that resummation schemes must be used. Moreover, we explain why such a model, 
based on the combination of perturbation theories and halo models, remains consistent with standard perturbation 
theory up to the order of the resummation scheme. We describe a simple implementation, based on a 1-loop "direct 
steepest-descent" resummation for the 2-halo contribution that allows fast numerical computations, and we check that 
we obtain a good match to simulations at low and high k. We also study the dependence of such predictions on the 
details of the underlying model, such as the choice of the perturbative resummation scheme or the properties of halo 
profiles. Our simple implementation already fares better than standard 1-loop perturbation theory on large scales and 
simple fits to the power spectrum at high k, with a typical accuracy of 1% on large scales and 10% on small scales. 
We obtain similar results for the two-point correlation function. However, there remains room for improvement on the 
transition scale between the 2-halo and 1-halo contributions, which may be the most difficult regime to describe. 

Key words, gravitation; cosmology: theory - large-scale structure of Universe 



1. Introduction 



In the standard cosmological scenario, the large-scale struc- 
tures we observe in the recent Universe (galaxies, clusters, 
filaments, voids, ..) have formed through the amplification 
by gravitational instability of small primordial perturba- 
tions genera ted in the pr imordial Universe by quantum 
fluctuations (|Peebleslll980() . Moreover, at the beginning of 
matter domination (z ~ 3000), the power inc reases on sma ll 
scales for cold dark matter (CDM) models (|Peebleslll982T) . 
which leads to a "hierarchical scenario" where smaller scales 
become nonlinear first. Then, increasingly large and mas- 
sive structures form in the course of time, as small struc- 
tures merge and larger scales turn nonlinear. Therefore, on 
large scales or at early times, it is possible to use linear the- 
ory or, more generally, perturbation theory, while on small 
scales or at late times one must use numerical simulations 
or phenomenological models. A quantity of great interest 
to characterize the density fields built by these processes 
is the density power spectrum P(k), which is the Fourier 
transform of the density two-point correlation function. On 
very large linear scales this is actually sufficient to fully de- 
termine the statistical properties of the matter distribution 
if the initial conditions are (almost) Gaussian, as in stan- 
dard scenarios of single-field inflation. On smaller nonlinear 
scales, higher order correlations are generated by the non- 



linearities of the dynamics, but P(k) remains a standard 
tool to compare observations with theoretical predictions 
(since higher order statistics are increasingly noisy or diffi- 
cult to predict). 

When one focuses on large scales, several observational 
prob es, such as w e ak le nsing surveys (M assev et al.l 
120071 : iMunshi et all I2008D . or galaxy redshift sur- 
veys provide rich cosmological information. For ex- 
ample, t he signature of baryo n acou stic oscillations 
(BAO) (lEisenstein et a .1 Il998l I2005T). redshi ft -space 
distortions (|Hamiltonl Il992t iPerc val fc White! l200l . 
non-Gaussianities in the primordial pertur bations 



(|Liguori et al.l 120101: Desiacaues k, Seliak 20101) and 



mass of neutrinos (jSwanson et alT 20101 : ISaito eTalll20Toh 



are sensitive to the clustering on linear to weakly nonlinear 
scales, where the departures from the linear regime are 
already noticeable but still modest. In order to meet the 
accuracy of future observations a theoretical accuracy on 
the order of 1% is required, which is beyond the reach 
of simple phenomenological models or fits to numerical 
simulations (if one wishes to obtain robust predictions 
for a large range of initial conditions and cosmological 
parameters). However, this regime should be within the 
range of validity of perturbation theories, which offer the 
advantage of systematic and reliable predictions, without 
fitting parameters. This has led to a renewed interest 
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in perturbative approaches in recent years, as it may 
be possi ble to improve ov e r the standard perturba tion 
theory (|Goroff et all 119861: iBernardeau et al.l l2002h by 
using resummation schemes that allow partial resum- 
mation s of higher order t erms ([Crocce fc Scoccimarrol 



2006 bJlal: IValaeeasI 2007at iMatarrese fc Pictroni 2007; 



Taruva fc Hiramatsul 120081: IValageasI 120081; iPietronil 120081: 
MatsubaraT 20081 ) . In particular, this somewhat improves 



the accuracy of the predicted matter power spectrum 
on the large scales associated with BAO, as compared 
with st andard perturbation theory t r uncated at the same 
order (ICrocce fc Scoccimarrol 120081: ICarlson et all 120091: 
iTaruvaet al.ll2009ft . 



On the other hand, many observational probes, such 
as weak lensing surveys on smaller angular scales and 
galaxy surveys, are sensitive to highly nonlinear scales. 
There, systematic analytical approaches that have been 
proposed so far do not apply, in particular because of 
the importance of shell crossing effects that are not ad- 
equately taken i nto account (an exception is the formal- 
ism developed in[Valagcas (2004), which does not make the 
single-stream approximation since it deals with the phase- 
space distribution /(x, v;t), but leads to heavy compu- 
tations). Then, one must resort to numerical simulations 
or to phenomenological models. One such model is the 
Lag rangian mapping p resented in lHamilton et al.l (|l99ll ) 
(see iPeacock fc Doddsl (|l996l ) in Fourier space) , that writes 
the nonlinear two-point correlation or power spectrum in 
terms of its linear counterpart on a different scale. A sec- 
ond model is the "halo mo del" (jScherrer fc Bertschingerl 
119911: ICoorav fc Shethl I2002D . where the matter distribu- 
tion is described as a collection of halos. This yields two 
contributions for the power spectrum or two-point correla- 
tion function, a "1-halo" term associated with particle pairs 
that belong to the same halo, and a "2-halo" term associ- 
ated with particles that belong to two different halos. Then, 
the 1-halo term, which dominates on small scales, is sensi- 
tive to the density profile and mass distribution of the halos, 
whereas the 2-halo term, which dominates on large scales, 
is sensitive to the correlation between halos. In practice, 
one often replaces the 2-halo term by the linear two-point 
correlation, or power spectrum, to reproduce linear the- 
ory on large scales . Prev i ous works (|Coorav fc Shethi r2002: 
ISmith et all 120031 l2007t iGiocoli et all I2010D haveshown 
that halo models provide a good match to numerical simu- 
lations, especially at high fc, whereas Lagrangian m appings 
based on the stable clustering ansatz (jPeeblesll 19801 ) do not 
fare as well. 

It is clear that it is of great interest to develop unified 
models, which combine the benefits of systematic pertur- 
bation theories on large scales with the reasonable match 
to simulations on small scales provided by halo models. 
The goal of the present paper is precisely to investigate 
the building of such a model, that describes all scales, from 
the linear to the highly nonlinear regime. We can already 
note here two points that must be addressed to reach this 
goal: i) the increasingly large growth at high k of higher 
order terms obtained within standard perturbation theory, 
and ii) the nonzero limit for k — > of the 1-halo term (as 
written in previous works), which eventually becomes non- 
negligible (and even dominant) on large scales, since the 
linear power spectrum typically scales as fi(fc) ~ fc at low 
fc for CDM scenarios. Thus, the standard implementations 



of both frameworks are badly behaved beyond their regime 
of validity, where they should actually become negligible. 

In this work we address these issues, and we perform a 
detailed comparison between different possible implemen- 
tations of such a unified model, as well as with numerical 
simulations, as follows. We first consider the computation 
of the density power spectrum from a Lagrangian point of 
view in Sect. O This provides a decomposition over 2-halo 
and 1-halo terms that are slightly different from the versions 
found in previous works, and we emphasize their relation- 
ship with "perturbative" and "non-perturbative" terms. In 
particular, we explain why the 1-halo contribution contains 
a counterterm, missed in previous studies, that ensures that 
it becomes truly negligible on large scales, solving the sec- 
ond point ii) above. Next, we present a simple implementa- 
tion of this model in Sect. [31 and we explain how the point 
i) above is solved by using resummation schemes instead of 
the standard perturbation theory. We describe in Sect. [4] 
the N-body simulations that we use to evaluate the accu- 
racy of our models, and we perform detailed comparisons 
for the nonlinear density power spectrum in Sect. [3 from 
k ~ 0.02/iMpc -1 up to k ~ 100/iMpc" 1 , and for redshifts 
z = 0.35 up to z — 3. Then, we investigate in Sect.[B]the de- 
pendence of our results on various ingredients of the model, 
such as the properties of virialized halos or the choice of the 
perturbative resummation scheme. We also evaluate the im- 
portance of the 1-halo counterterm on large scales. Next, 
we consider the real-space two-point correlation function 
in Sect. [7] Finally, we estimate in Sect. [5] the accuracy of 
such unified models, for both the power spectrum and the 
two-point correlation, and we conclude in Sect. [5J 

In this paper we neglect the impact of baryonic physics 
on the matter power spectrum, or more precisely we do not 
explicitly address this issue. Thus, we assume that it can 
be incorporated in an effective manner through the choice 
of the halo density profile that enters the model. 



2. Decomposition of the density power spectrum 
from a Lagrangian point of view 

We show in this section how the density power spectrum 
can be split into "perturbative" and "non-perturbative" 
terms, and how they are related to the 2-halo and 1-halo 
terms, from a Lagrangian point of view. 

2.1. "Perturbative" and "non-perturbative" terms 

Let us recall that in a Lagrangian framework one considers 
the trajectories x(q, t) of all particles, of initial Lagrangian 
coordinates q and Eulcrian coordinates x at time t. In par- 
ticular, at any given time t, this defines a mapping, q n> x, 
from Lagrangian to Eulerian space, which fully determines 
the Eulerian density field p(x) through the conservation of 
matter, 



/j(x) dx = pdq, 



(1) 



where p is the mean comoving matter density of the 
Universe and we work in comoving coordinates. Then, 
defining the density contrast as 



p{x,t) 



(2) 
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and its Fourier transform as 
dx 



,-ik-: 



^ J W 
we obtain from Eq.([T]) 



5(x), 



~S(k) 



_dq_ 

(2tt) 3 



e -ik-x(q) _ g-ik q 



(3) 



(4) 



this stan dard perturbative series have b een recently in- 
troduced (ICrocce fe Scoccimarrol 2006bllal: IValaee as 2007a; 



Next, defining the density power spectrum as 
(*(ki)S(k 2 )> = S D (ki + k 2 )P(fc 1 ), (5) 
we obtain from Eq.Q, using statistical homogeneity, 



P(k) 



(0) 



where we introduced the Eulerian-space separation Ax, 

Ax = x(q) - x(0). (7) 

The expression (j6]) is fully general since it is a simple con- 
sequence of the matter conservation Eq.([T]), and of statis- 
tical homogeneity. In particular, it holds for any dynamics, 
such as the one asso ciated with the Zeldovich approxima- 
tion (Zcldovichl l 19701) . There, the mapping x(q) is given by 
the linear displacement field as x(q) = q + \J/i(q), so 
that for Gaussian initial conditions (i.e. >&L(q) is Gaussian) 
one can easily compute the average ^ and derive an ex- 
plicit analytic expression for th e density power spe ctrum 
([Schneider fe Bartelmaiml Il995t iTavlor fe Hamilton! [19961 
!Valageasll2010a|) . 

The second term in Eq.©, e lkq , merely gives a Dirac 
factor, <5n(k), and can often be discarded in exact or sys- 
tematic computations (e.g., perturbative expansions) if one 
consider k ^ 0. Within perturbative expansions of P(k) 
over powers of the linear power spectrum, P^(k), it simply 
cancels out the zeroth-order term so that P(k) = PL(k) at 
lowest order. However, as we shall see below, it will prove 
important in the following analysis to keep track of this 
factor. Indeed, it will be split over two parts, along with 
the first term e lkAx , and it is clear that one should either 
discard it or keep it in a consistent manner in both parts, 
contrary to what is implicitly done in usual versions of the 
halo model. 

For the three-dimensional gravitational dynamics it is 
not possible to derive an explicit analytical expression for 
the average Therefore, as in the Euleri an framework, 
one must resort to pertu rbative expansions (jBouchet et al.l 
Il995t iMatsubaral I2008D or phenomenological models. In 
this paper, our goal is precisely to use Eq.® as a start- 
ing point to decompose the density power spectrum, P{k) 1 
over two parts, associated with "perturbative" and "non- 
perturbative" contributions. 

Here and in the following, we call "perturbative" 
those contributions that can be obtained from pertur- 
bative expansions over powers of the linear power spec- 
trum Ph(k) within the fluid approximation. In practice, 
these may be derived from the equations of motion, writ- 
ten in terms of the density and velocity fields in the 
usual hydrodynamical description of the system, by writ- 
ing these nonlinear Eulerian fields as pert urbative se- 
ries over powers of the l i near growing mode (jGoroff et al.l 
1986; iBernardeau et al.l [2002). Partial resummations of 



Matarrcs efe Pietronil 120071: iTaruva fe Hiramatsul 12008 
IValaeeasI 120081: iPietronil I2008D but remain within this hy 
drodynamical regime. Of course, it is possible to develop 
perturbative expa nsions, and associated resummations, i n 
Lagrangian space (jBouchet et al.| [T995: Matsubara 2008), 
by looking for an expansion of the displacement field ^(q) 
over powers of the linear displacement field. Taken at face 
value, these schemes imply some shell crossing. For in- 
stance, at lowest order one recovers the Zeldovich dynamics 
for the displacement field. However, this does not mean 
that they take into account in a systematic fashion the 
physics beyond shell crossing, which remains beyond the 
reach of current approaches of this kind. Technically, the 
gravitational potential is obtained from the Poisson equa- 
tion, where the density contrast is obtained from the conser- 
vation of matter ([T]) as <5(x) = | dct(3x/9q)| _1 — 1, and one 
expands the Jacobian over powers of the displacement 
In doing so, one discards the absolute value and misses the 
changes of sign that occur after shell crossing. On the other 
hand, these schemes give the same final expansion over pow- 
ers of Pl for the nonlinear power spectrum as the Eulerian 
schemes, up to the order of truncation of the computa- 
tion. Therefore, we consider both perturbative approaches 
as "perturbative" schemes, in the sense of expansions over 
powers of Pl that neglect shell crossing. By contrast, we 
call "non-perturbative" those contributions that cannot be 
obtained through these approaches, either because they are 
beyond the reach of perturbative expansions (such as fac- 
tors of the form e~ 1 l Ph that cannot be uncovered by Taylor 
series over Pl) or because they arise from the behavior of 
the system beyond shell crossing (and require going beyond 
the single-stream approximation). 



In IValaeeasI (|2010al) such a decomposition between per- 
turbative and non-perturbative contributions to the den- 
sity power spectrum was obtained within the framework of 
the Zeldovich dynamics. More precisely, a "sticky model" 
that only differs from the Zeldovich dynamics beyond shell 
crossing was introduced. Then, whereas the Zeldovich den- 
sity power spectrum is exactly given by a resummation of 
perturbative terms (i.e., the non-perturbative contribution 
is zero), the "sticky model" power spectrum contains an 
additional non-perturbative contribution, while its pertur- 
bative part is identical to the one of the Zeldovich dy- 
namics. Inspired by this simpler example, we perform a 
similar decomposition of Eq.© over perturbative and non- 
perturbative contributions. However, we no longer have a 
simple criterion to distinguish particle pairs that have un- 
dergone shell crossing, because for the gravitational dynam- 
ics even at the perturbative level (at third order) the dis- 
placement field develops rota tional terms (|Buchert 1994; 
IBernardeau fc Valaeeasl [20081 ). Therefore, we use a phe- 
nomenological halo model, in the sense that we assume that 
the matter distribution can be described as a collection of 
spherical halos defined by a given nonlinear density contrast 
6. As usual, we can write the halo mass function as 



n(M)dM =■£/(!/) — , withi/=-^ ; 

M v <t{M ) 



(8) 
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in terms of a scaling function f(v) and of the reduced vari- 
able v. Here, a(M) is the rms linear density contrast at 
scale M, or Lagrangian radius qM, with 



47T 

cr(M) = a(q M ) with M = p—q M , 
and 



a\q) = AK dkk^P L {k)W{kqY 
Jo 



(9) 



(10) 



where W(kq) is the Fourier transform of the top-hat of 
radius q, defined as 

W(kq) = [ £ = 3 Sin(fcg) ~"f° s(fcg) . (11) 

In the second Eq. ([8]) the linear density contrast 6l is related 
to the nonlinear density threshold 5 that defines the halo s 
through the spherical collapse dynamics (Valageas 2009a), 



(12) 



and for Gaussian initial conditions the large-mass tail shows 
the exponential falloff 



v -> oo : ln[/(z/)] 

whence 

M oo : n(M) ~ er 6l/ 



(13) 
(14) 



This corresponds to the Press-Schechter falloff 
(jPress fc Schechterl Il974tl . but with a somewhat lower 
threshold 8l given by Eq. p^|) (the usual Press-Schechter 
threshold actually corresponds to Sl = .F _1 (oo), associ- 
ated with full collapse to a point). However, the relation 
(fT2")) only holds for S < S vir , which typically gives 5 < 200, 
as larger density contrasts are associated with inner 
shells where shell crossing plays a key role and modifies 
Eq.dTz ]!). associated w ith spherical dynamics at constant 
mass (Valageas 2009a). On the other hand, the halo mass 
function satisfies the normalization 



poo j 



(15) 



which ensures that all the mass is contained within such 
halos (for linear power spectra such that a(q) grows to in- 
finity on small scales). This also ensures that there is no 
overcounting (as would be the case if one used a mass func- 
tion with a normalization greater than unity). 

Then, in Lagrangian space, the probability dF for one 
particle qi to belong to a halo of mass in the range [M, M+ 
AM] reads as 



(16) 



This is also the fraction of matter enclosed within such 
halos. Next, making the approximation that each halo 
comes from an initial spherical region in Lagrangian space, 
the probability for a second particle q2, at distance q = 
1 12 — qil, to belong to the same halo reads as 



Fm ^> = 4, 3 7Z 



V) 



(17) 



Here we integrate over all positions qi within the spherical 
volume V of radius qM, and we integrate over all directions 
Vl q of the Lagrangian vector q = q2 — qi- The top-hat 
factor 0(q2 6 V) is unity if q2 is located within the volume 
V, and zero otherwise. By isotropy the result only depends 
on the length q, and performing the integrations yields 



< q < 2q 



M 



F M (q) 



(2q M -g) 2 (4g M + q) 



(18) 



and Fm(q) = for q > 2qM- Therefore, combining Eas. (Tl7o| 
and ([18]) . we obtain the probability that a pair of particles 
of Lagrangian separation q belongs to the same halo of mass 
M as 



q M >-: dF q (M) = _ /(,)-, (19) 



and dF q (M) = for q^j < q/2. In particular, the proba- 
bility that the pair {qi,q2} belongs to a single halo writes 
as 



Fm(q) 



°° (2q M ~q) 2 (MM + q) 



V q /2 



(20) 



where v q /2 is defined as in Eq.®, for the Lagrangian ra- 
dius qM = q/2. On the other hand, the probability F 2 ]i{q) 
that the Lagrangian pair does not belong to a single halo 
(whence the two points belong to two different halos) reads 
as 



F 2 H{q) = 1-Fm(q). 



(21) 



Then, we split the average in Eq.® over two terms, 
Pih and P2H, associated with pairs {0, q} that belong to a 
single halo or to two different halos, 



P(k) = P m (k) + Fm(k), 
with 



Pm(k) = J -0^F 1R (q) (e ik - Ax -e ik "») 



1H 



and 



P2u{k) = J J^FMq) (e ik Ax -e ik q ) 2 H. 



(22) 



(23) 



(24) 



Here the averages (...)ih and (...)2H are the conditional av- 
erages, knowing that the pair of length q belongs to a single 
halo or to two halos. The decomposition (p?!?|) clearly cor- 
respon ds to the 1-halo and 2- halo terms of the usual halo 
model (|Coorav fc Shethll2002t) . Then, to make the connec- 
tion with the distinction between perturbative and non- 
perturbative terms, we note that at a perturbative level 
-FiH is identically zero and P2H unity, 

Fm = and F 2 h = 1, 

at all orders of perturbation theory. (25) 

Indeed, the large-mass tail (fl4")) is actually a rare-event limit 
that holds both in the large- mass limit , at fixe d linear den- 
sity power spectrum J Valageas' l2002bl . l2009al) . and in the 
quasi-linear limit at fixed mass, where the amplitude of the 
linear density power spectrum goes to zero (|Valageasl2002al 
l2009al ). It is this second regime which corresponds to usual 
perturbation theories, where as recalled above we look for 
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expansions over powers of the amplitude of the linear power 
spectrum. Then, because of the exponential decay of the 
form e~ x l a ( M > we can see that the expansion over powers 
of Pl of Fm(q) defined in Eg . (|2"0)) . at fixed q, is identically 
zero. From Eq. (l2"TT) this also yields F 2 h = 1 at all orders of 
perturbation theory. Therefore, we can see from Eqs. (|23[) - 
that the 1-halo contribution is a fully non-perturbative 
term, while the 2-halo contribution is (almost) the pertur- 
bative term multiplied by the factor F 2 n(q). 

The factor Fm being non-perturbative is not mainly 
related to shell crossing, but simply to the fact that it can- 
not be recovered by a series expansion over powers of Pl ■ 
However, for halos that would be defined by a high density 
threshold, typically 8 > 200, the exponential falloff (TH)) 
is modified by shell crossing ( i .e. the factor 8l is no longer 
given by Eq. (IT2l . see Valageas (20093)) so that it would also 
be non-perturbative in this sense. On the other hand, even 
for lower threshold 8, the exact form of the factor Fm(q), 
and in particular the low-mass tail of the halo mass func- 
tion, depends on the behavior of the system beyond shell 
crossing. 

The decomposition (|2"2"]l . with the Lagrangian-based in- 
terpretation (l23l) - (j24)) . has the advantage to automatically 
satisfy the conservation of matter. Thus, thanks to Eq. (|2"Tj) 
we count all particle pairs once. By contrast, in the Eulerian 
derivation of the halo model, where we first write the den- 
sity field jo(x) as a sum of halo profiles, we would need 
to pay attention to possible overlaps between halos, which 
arise when we use a spherical approximation. Thus, the 
splitting (|22p is more easily expressed in this framework, 
and one can independently focus on the modelization of 
the averages (e lkAx — e lk q ). The latter also offers a closer 
link to the dynamics, through the mapping q M> x. We shall 
not make much use of this relationship in the following, as 
we use simple approximations that allow us to recover the 
usual Eulerian expressions, with the addition of simple pref- 
actors and counterterms, but this may provide a route to 
more accurate modeling. 

2.2. "2-halo" contribution 

We first consider the 2-halo contribution (|24[) . Since the 
conditional average involves the constraint that the pair 
{0, q} does not belong to the same halo, the mean of e lk ' Ax 
is "biased" as compared with a mean over all possible pairs. 
However, in order to simplify the computation of this term, 
we note that at a perturbative level F 2 r = 1, as seen in (1231) . 
so that we can replace the average (...)2H by the mean over 
all pairs, as given by perturbation theory, 

P 2R (k) ~ J ^3 F 2H (g) (e ik Ax - e ik <% ort , (26) 

where the subscript "pert" denotes quantities obtained 
from standard perturbation theory (or equivalently its vari- 
ous resummation schemes) . Thus Eq. (l26l) is still exact at all 
orders of perturbation theory. This expression is best suited 
for perturbation theories developed within the Lagrangian 
framework. Unfortunately, as we shall discuss in Sect. 16.51 
below, Lagrangian perturbation theories built so far are not 
as efficient as their Eulerian counterparts, especially when 
we consider available resummation schemes. Then, in or- 
der to make contact with Eulerian perturbation theories 
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we further approximate Eq. (|26[) as 

P 2H (fc) * F m (l/k) J (0!<e ik ' Ax - expert 

= F 2H (l/k)P pcrt (k), (27) 

where we have replaced the q— dependent factor F2n(q) by 
its value at a typical scale q ~ 1/k. Again, this is legitimate 
at a perturbative level, where F 2 h = 1, so that Eq. (|2"?| re- 
mains exact at all orders of perturbation theory. To obtain 
the second line we simply used the exact expression ©, 
which implies the same equality in terms of perturbative 
expansions. 

Let us recall that from (|2"5"|) the 1-halo contribution is 
zero at all orders of perturbation theory, so that the full 
power spectrum P(k) of Eo. (j22"|) automatically agrees with 
perturbation theory at all orders, whether we use Eq. ([23| 
or Eq. (|27p . Then, within these approximations the only ef- 
fect of non-perturbative corrections to the 2-halo term is to 
multiply the perturbative power spectrum P per t(fc) by the 
prefactor F 2 n(l/k). 

Within the usual h alo model the 2-halo term reads as 
(|Coorav fc Shethll2002ft 

Pm m (k) = I ^/(fi)/h)% 1 (l:)iiM 2 (fc)PM 1 M 2 (fc) 
~ P L (k), (28) 

where Uj\t{k) is the normalized halo density profile, defined 
in Eq. (pn"j) below, and PM 1 M 2 {k) is the halo power spec- 
trum. The second line (|28[) is obtained in the low-fc limit, 
so that u M (k) -> 1 and P Ml M 2 (k) ~ b{M x )b{M 2 )P L {k) , 
with a halo bias b(M) that is normalized to unity. It is 
also possible to combine perturbation theory and nonlin- 
ear halo bias to make the expression above consistent with 
standard perturbation theory whi le building a mode l for 
the halo power spectrum itself, see lSmith et al.l (|2007h . 

Of course, in order to describe the weakly nonlinear 
regime one can as well replace -Pi(fc) by P per t(k) in Eq. (|2"8"|) . 
which gives an expression very similar to Eq. (l27l) . Then, 
we can see that a first difference between the halo-model 
expression (|28fl and Eq. (l2~Tl) is that we did not need to in- 
troduce any halo bias to derive Eq. (|2"T|) . In fact, although 
the contribution ([24| is associated with a 2-halo term, by 
making the simple approximation (|26p we avoid any need 
to consider halo biasing, and as explained above this does 
not spoil the agreement with perturbation theory. The sec- 
ond difference is the prefactor F 2 n{l/k) in Eq. (l2~T|) . As 
seen from Eq.© and the splitting (|2"2"|) . this term is re- 
quired by self-consistency, to ensure that the two averages 
*1H(«) (e ik Ax )m and F m (q) (e ik Ax ) 2H sum up to (e ik Ax ), 
and in particular that they sum up to unity for k — > 0. 
Within the usual halo model, this factor is implicitly set 
to unity by taking the large-scale limit in Eq. (j2"81 and ig- 
noring exclusion constraints on the halos. This is also valid 
within perturbation theory, as seen in (125|) . However, in 
order to describe the weakly nonlinear regime, where the 
1-halo term is nonzero, it is best to keep the prefactor 
F 2 n(l/k) in Eo. ([27|) . to keep a consistent model and to 
avoid any overcounting. 



6 



P. Valageas and T. Nishimichi: Combining perturbation theories with halo models 



2.3. "1-halo" contribution 



From Eas. (j20[) and (123)) the 1-halo contribution to the den- 
sity power spectrum reads as 



f dq f°° d» 
Pm(k) = / — ^ / — / 



(2 gM -g) 2 (4g M + g) 

16^ 



x(e i 



k-Ax 



,ikq 



(29) 



In order to compute the average in Eq. (|29[) . within a halo 
of mass M, we make the approximation of fully virialized 
spherical halos. Thus, we describe the halos as spherical ob- 
jects, truncated at a radius tm such that the mean density 
contrast within this radius is the nonlinear threshold S used 
to define these objects, and the enclosed mass is equal to 
M, 



(30) 



47T o 4-7T o 

M = p —q M = (1 + S)p—r 3 M . 

We introduce as usual the normalized Fourier transform of 
this halo radial profile, 



uu(k) = 



Jdxe iU ' x pAi(x) _ 1 



/ dxp M (a 



M 



p M (x), (31) 



where pm{x) is the halo density profile. Next, using the 
approximation of fully virialized halos, that is, that the 
two Lagrangian particles "0" and "q" have lost all memory 
of their initial locations and are independently located at 
random within the halo, we write 



/ ik-Ax\ 



l - J d Xl dx 2 PM (xiVA f (x 2 ) e ik -( X2 - Xl >(32) 

(33) 



M 2 _ 
u M {k) 2 



Substituting into Eq. (|2"9"]) and exchanging the order of in- 
tegration gives 



r°° dv r 2qM 

Pm(k) = / —f{v) / 

Jo v JO 



dq {2q M - q) 2 (MM + q) 



(2tt) 3 



16qlj 



x (u M (fc) 2 -e ik -i), 
and the integration over q yields 

poo j AT/ 

Pm(k) = / —f{v)——- (u M (kf - W(kq M f ) . (35) 



(34) 



o ^ p(2tt)= 

Theref ore, we recover the us ual 1-halo term of the halo 
model (jCoorav fc Shethll2002tl , with the addition of the new 
counterterm W(kqM) 2 , where W was defined in Eq . (fTTl) . 

This counterterm originates from the second term in 
Eq.([6]), which actually sums up to a Dirac factor #o(k). 
However, since we eventually split the nonlinear power 
spectrum as in Eq. (|2"2")l and we use different approxima- 
tions for the 2-halo and 1-halo contributions, it is best to 
keep track of this factor in both contributions. In particu- 
lar, as pointed out in Sect. 12.11 one should avoid keeping 
it in one contribution and disregarding it in the other one. 
Thus, it has been explicitly used in Eq. (|2"T|) to cancel out 
the zeroth-order term F 2 H<!>D(k), so that one should keep 
it in the 1-halo term f|23|) . Within the usual halo model 
this term is usually missed because the splitting between 
the 2-halo and 1-halo terms is not treated so carefully. In 
particular, as noticed in Sect. 12.21 the usual halo model 
implicitly takes Fzn = 1 (which is valid in perturbation 



theory) while taking Fm ^ by adding the 1-halo term, 
which is not fully consistent. 

In fact, the counterterm of Eq. (j35|) can be recovered 
without going through the previous steps, associated with 
a Lagrangian point of view, by the following simple ar- 
gument. Let us consider a uniform universe, where there 
are no density fluctuations and the density is everywhere 
equal to p- Then, within the spirit of the usual (Eulerian) 
halo model, and making the same approximations (e.g., ne- 
glecting geometrical constraints, associated with exclusion 
constraints, and departures from spherical symmetry), we 
can consider that all the matter is within halos of constant 
density p with an arbitrary distribution of halo radii. That 
is, we can split this uniform system over an arbitrary set 
of cells, which we can call halos. Then, within such a halo 
model the usual 1-halo term for the full density correlation, 
Pgg(fc) (defined as (p(k x )p(k 2 )) = tf D (ki + k 2 )P""(fci), as 
in Eq.©), reads as the counterterm of Eo. ((35|) . with a posi- 
tive sign and multiplied by p 2 . Indeed, the density profile of 
such halos is simply the top- hat of radius qM, since 5 = 
whence ru = qM, so that the normalized Fourier trans- 
form of the halo radial profile defined in Eq. (|3"T1) is equal to 
the Fourier transform (fTTj) . Ujvr(fc) == W(kqnj)- This means 
that within a halo model, for any halo mass function, the 

1- halo contribution to the power spectrum associated with 
an uniform-density medium reads as the counterterm of 
Eq. (|35"]) . Since the density-contrast power spectrum can also 
be defined as p 2 P{k) = P^(fc) - P^fc) (i.e. we consider 
the connected two-point correlation), we must subtract this 
uniform-density term as in Eq. (|35[) . as the first term UM(k) 2 
corresponds to the total density power spectrum, P pp (k). 

In agreement with this result, we can note by an in- 
tegration over k of the second term of Eq. (f3~i)) that the 
second term of Eq. (|35[) is indeed a well-normalized approx- 
imation to 5o(k), where the Dirac peak is broadened over a 
size ~ 1/qM associated with the typical halo scale (thus, as 
could be expected, this becomes an increasingly good ap- 
proximation to <5o(k) as the halo size goes to infinity and 
the associated power is repelled to k —> 0). Clearly, the fact 
that this does not yield exactly a Dirac factor, whence that 
it does not vanish for nonzero k, is due to the fact that it is 
only a partial contribution to the full power spectrum, since 
by construction within such halo models there is always a 

2- halo term, and it is only the sum of both contributions 
that contains such a Dirac factor. (This also agrees with the 
fact that we recover a Dirac factor in the limit qM — > oo, 
since at fixed scale 1/k the 2-halo term clearly disappears 
in this limit). 

In practice, the counterterm of Eq. (|3"5f is small on 
most scales of interest, especially in the nonlinear regime. 
However, it plays an important role on very large scales. 
Indeed, since «m(0) = 1 we can see that the first term 
of Eo. ((35|) . associated with the usual halo model, goes to a 
nonzero constant as k — > 0. Since for CDM-like initial power 
spectra Ph{k) oc k Us with n s ~ 1 as k — >• 0, this implies 
that the usual 1-halo term dominates on very large scale, 
which is not correct. In fact, from the conservation of mat- 
ter a small-scale redistribution of matter generically yields 
a P(k) oc k 2 tail at low k, while taking int o account mo - 
mentum conservation one obtains a fc 4 tail (|Peebleslll974[ ). 
To solve this problem the use of compens ated halo profiles 
(i.e. w ith «m(0) = 0) was investigated in ICoorav k, Shetbl 
(|2002D . However, they noticed that this recipe fails because 
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it also spoils the 2-halo term, as seen from the first line in 
Ea. (|2"5|) . Here we can see that this problem is automatically 
solved within our approach by the counterterm W(kqM) 2 
in Eg. (1351) . In particular we recover at low k the expected 
slope Pm(k) oc fc 2 , associated with small-scale rearrange- 
ments (indeed, the I-halo term only describes the redis- 
tribution within small clumps, and is not sensitive to the 
long-range correlations between clumps, that is described 
by the 2-halo term). 

Looking at the steps of the derivation of Ea. (f35| we 
can also see why we fail to recover the fc tail due to mo- 
mentum conservation. Indeed, in the approximation (1331) . 
where we have assumed full virialization, we have erased 
all memory of the initial positions and velocities of the two 
Lagrangian particles {0, q}, which clearly implies that we 
have disregarded the constraints associated with momen- 
tum conservation. Therefore, in order to obtain a fc 4 tail 
one should improve this approximation, by taking care at 
some level of the particle velocities. Of course, a simpler 
remedy would be to modify the counterterm W(kqMj in 
Ea. (|3"5"|) . by using a window H^fc^A/) 2 that cancels both 
the terms k° and k 2 of mm (fc) 2 - However, since the CDM 
linear power spectrum roughly behaves as Pl (fc) oc fc at low 
fc the fc 2 tail is sufficient to make the 1-halo term subdom- 
inant on large scales, as we shall check in Sect. \5\ below. 
Therefore, we shall keep the simple approximation (|33p in 
the following, and Ea. (|3"5"|) for the 1-halo term, as it appears 
to be sufficient to reach a good agreement with numerical 
simulations. 



2.4. Extension to the galaxy power spectrum 

Although in this article we focus on the dark matter power 
spectrum, we briefly discuss in this section how the previous 
results can be used for the galaxy power spectrum, which 
is also a quantity of great practical interest. This necessar- 
ily involves additional ingredients, as we must define the 
relationship betwe en the dark matter density field and the 
galaxies. As usual fSe lia3l2000t ICoorav fc Shethll2002t) . we 
keep the splitting (f22j) over 1-halo and 2-halo terms. Then, 
we simply write the galaxy 2-halo contribution as 



Pf^fc) = (b) 2 P m (k), 



(36) 



where (6) is the mean bias of the galaxy population that we 
consider. This can be computed using one of t he bias models 
that have been proposed in p revious work s (iMo fc W hite 
19961 ISheth fc Tormenl [llM iSheth etafl 120011: iValaeeas 
2009aL l2010rJ or fits to numerical simulations (Ti nker et al.l 
2010). Next, following the Eulerian interpretation of the 
counterterm of Eg. (|35[1 discussed in the previous section, 
as arising from the difference (n gal n gal ) — (n gal )(n gal ), we 
write the 1-halo term as 



<(fc) 



M -f_ (N(N-l)) 



M 2 



X U M (k) p -W(kq M )P 



(37) 



where n is the mean galaxy number density, (N(N — 1)) 
is the mean squared number of galaxies (minus a factor 
N) in halos of mass M, and p = 1 or 2 depending on 
whether the f ormer quantity is small er or greater than unity 
(|Seriakll2000HCoorav fc Shethll2002D . Thus, we have simply 



subtracted to the usual expression the counterterm asso- 
ciated with galaxies set at random within the Lagrangian 
radius qM (i.e., corresponding to a uniform density uni- 
verse). Again, the quantity {N(N — 1)) must be computed 
from a model defined for the population of galaxies that we 
consider. 

The derivation of Eg. (1371) is more phenomenological 
than the derivation of Ea. (|3"5"|) from the Lagrangian point 
of view (1291) . However, since one always needs to introduce 
some phenomenological model to relate galaxies to the mat- 
ter distribution (here through the two quantities (b)(M) 
and (N(N - 1))(M)), this should be sufficient. 

The contribution (|37l) to the galaxy power spectrum 
vanishes as fc 2 in the large-scale limit, which again solves 
the unphysical nonzero limit obtained in previous imple- 
mentations. This behavior remains valid even though galax- 
ies are discrete objects, since this property only implies 
a constant shot- noise asymptote for P gal (fc) in the high- 
fc limit. This shot-noise has actually been subtracted in 
Eq.fl2Z}, using the quantity ( N(N - 1)) instead of (AT 2 ), 
to follow the usual convention (|Peeblesj|1980l; ISeliakll2000t 
ICoorav fc Sheth|[200l . At low fc, we are dominated by the 
2-halo contribution and in case of constant large-scale bias 
we recover the slope P gal (fc) oc fc™ s of the primordial matter 
power spectrumQ 



3. A simple implementation 

We describe in more detail in this section how we implement 
Eos. (f2"7f and ((35)) to compute the nonlinear density power 
spectrum. 



3.1. "2-halo" contribution 

We first consider the 2-halo contribution (|27|). and more 
precisely the perturbative part P pert (fc), as we shall discuss 
the non-perturbative prefactor F2n(l/h) in Sect. l3.2l below. 

Since we cannot compute and resum all terms of 
the perturbative expansion of P(k), contrary to the 
simpler Zeldovich dynamics, we must use for P pG rt(fc) 
an approximation associated with a truncation at a fi- 
nite order. This leads to several possible choices, since 
a priori we can use the standard pertu rbation theory 
(|Goroff et al.l ll986; Bcr nardeau et al.ll2002ft as well as any 
of the various resummatio n schemes that have been intr o- 
duced in the recent years (ICrocce fc Scoccimarrol 2006b. a; 
| Valageasll2007al: iMatarrese fc Pietronill2007l ; IValageasll2008l: 
iTaruva fc Hiramatsull2008tlPietronill2008l ). Indeed, all such 
methods give expressions for P per t(fc) that agree up to the 
order of the computation, and only differ by higher order 
terms. 

However, it turns out that standard perturbation theory 
cannot suit our purposes b ecause higher order term s grow 
increasingly fast at high fc (Be rnardeau et aDl2002th As we 



1 A simple illustration of a density field made of discrete 
objects that o beys such behaviors is provided by the ad- 
hesion model (|Gurbatov et alj 119891 ; IVer gassola e t al.l 119941 ; 
iBernardeau fc Valageas! l2010bD . This yields a power spectrum 
with a universal fc° slope at high fc, because of the pointlike 
density pea ks, while the slope at low fc d epends on the initial 
conditions ( Valageas & Bcrnardeau 20TFj), and for the simple 
cases n s = and n a = — 2 in ID one can obtain its exact ex- 
pression <|Valag cas 20091©. 
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shall check in Sect. 16.41 below, at one- loop order it already 
gives a steep contribution that remains non-negligible at 
high k, because the prefactor F2n(l/k) does not go to zero 
very fast (because for CDM power spectra a(q) only grows 
logarithmically at small q if n s = 1, and even remains finite 
if n s < 1). This implies that to get a well-behaved 2-halo 
term we need to use other schemes, that provide a perturba- 
tive part P pe rt(k) that remains small at high k as compared 
with the 1-halo term. Therefore, we must use one of the re- 
summation schemes that have been recently developed and 
show a well-behaved high-/c behavior. 

It is interesting to point out that, although such ap- 
proaches were devised in order to improve the accuracy of 
perturbation theory on large scale and to increase its range 
of validity (by improving its convergence through partial 
resummations of an infinite number of higher order terms) , 
our study shows that a second important benefit of these 
schemes is to provide well-behaved expressions for -P por t(fc). 
Here the point is not that the high-fc behavior is accurately 
reproduced (since anyway perturbation theory does not ap- 
ply in this range) but that it remains within reasonable 
bounds (typically close to or smaller than the linear power) 
so as to be subdominant. Thus, by giving systematic ex- 
pansions of Ppcstik), which agree with perturbation theory 
up to the order of the computation while remaining well- 
behaved at high k, these schemes allow the construction of 
unified models such as those studied in this article, based 
on Eq.(22), that show a correct behavior at both low and 
high k. (In the same spirit, we have checked in Sect. 12.31 
that the 1-halo term is well-behaved at low k.) This pro- 
vides a second important motivation for such resummation 
schemes, as compared with standard perturbation theory. 

In this article we fo cus on th e "direc t steepest-descent 
method" introduced in IValaeeasI (|2007al ). In terms of dia- 
grams it means that the two-point correlation C is given 
at one-loop order by the resummations shown in Fig. [TJ 
where single lines are the linear correlation and response 
functions Cl and Rl, while the double lines are the non- 
linear correlation and response f unctio ns C and R obtained 
at that order, see also lValagead (|2008l) . At this one-loop or- 
der the standard perturbation theory consists in keeping 
only the three diagrams with zero or one bubble among the 
infinite series shown in the second equality in Fig. [T] We 
recall the details of the computation of the (perturbative) 
density power spectrum -P per t(&) within this approach in 
App.El 

This direct steepest-descent method is not necessar- 
ily the most accurate resummation scheme. In particu- 
lar, it yields a response function that does not decay 
at high k or late times, but shows increasingly fast os- 
cillations with an amplitude that follows the linear re- 
sponse function. This is not realistic, since one expects 
a Gaussian-like decay for Eulerian response functions, as 
can be se en from theoretical arguments an d numerical sim- 
ulations (ICrocce fc Scoccimarrol I2006bllat iValageasI [20075 
iBernardeau fc Valageas! l2010al ). However, the fast oscilla- 
tions still provide an effective damping in a weak sense 
(that is when the response function is integrated over) . The 
reason why we consider this direct steepest-descent method 
here is that it provides a simple and efficient method, which 
satisfies our requirements and proves to be reasonably ac- 
curate, as we shall check in Sect. [5] below. Indeed, while by 
construction it agrees with standard perturbation theory at 
one-loop order, it is well-behaved at high k as it remains 



C = =<=C,=^= + 2 




♦-o-o<>o<> 

Fig. 1. The resummation performed by the "direct 
steepest-descent" method at one-loop order, for the two- 
point correlation C. The blue single lines are the linear 
correlation Cl and the red single lines are the linear re- 
sponse (propagator) Rl- The double red lines in the first 
equality are the nonlinear response R (at that order), which 
contains an infinite series of bubble diagrams, and explicit 
substitution gives the series of diagrams shown in the sec- 
ond equality below. 

close to the linear power spectrum on all scales (when we 
truncate the computation at one-loop order, as in this ar- 
ticle), see lValageasl (|2007al) . 

It is also particularly efficient because the integrals in- 
volved in the computation of the power spectrum factorize 
(because the bubble in the upper right diagram of Fig. [T] 
involves a product of linear two-point correlation functions, 
whose time-dependence can be factorized). This factoriza- 
tion can be clearly seen in Eq. (|A.29p . as the two-time cor- 
relation C(fc; 771,772) can be written as the sum of three 
products, each of them being of the form ^(771 ) x ^(772). 
Therefore, there is no need to compute two-dimensional 
time integrals, over both 771 and 772- A second property is 
that the integro-differential equations satisfied by the non- 
linear response i?(fc; 771, 772) (which is needed at intermedi- 
ate stages to obtain the power spectrum) can be reduced 
to (third-order) ordinary differential equations, as seen in 
Eqs. (|A.2ip - (|A.22p . This avoids the need to compute inte- 
grals over past times at each time step over 771 . These two 
properties allow a fast numerical computation, which can 
be a useful feature for practical purposes. 

3.2. "1-halo" contribution 

We now turn to the 1-halo contribution (|35D , which requires 
a model for the halo mass function and for the halo density 
profiles. In this article we consider the scaling function f(v) 
given by 

f{y) = 0.502 [(0.6 vf 5 + (0.62 is) - 5 ] e~ v °' /2 . (38) 

It satisfies the normalization constraint (1151) and it has al- 
ready been shown to provide a good match to numerical 
simulations for h alos defined by a density contrast <5 = 200 
(|Valageasll2009al) . Thus, substituting Eq.® into Eq.© 
we obtain the mass function of these halos, with a linear 
threshold S L = J r_1 (200) in the second Eq.©. At redshift 
z = this gives Sl — 1-59 for a ACDM cos mology, and it in- 
creases slightly to Sl — 1-6 at high z, see IValageasI (|2009al) 
and the fit of Eq.(12) therein. Next, the halo mass function 
also gives the probabilities Fm(q) and F2n(q) defined in 
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Eas.(j2"0 j) -(|2"T ]) . This yields in turn the prefactor F 2 a(l/k) 
that enters the 2-halo term (j2"Tj) . 

For the halo density profile we use the usual NFW pro- 
file (|Navarro et alj|1997| ). 



Pm{x) 



x/r s (l + x/r s y 



(39) 



which we truncate at the radius ^200, associated with the 
overall density contrast S = 200. The scale radius r s is 
defined in terms of the outer radius r2oo through the con- 
centration c(Af2oo), 



c(M 200 ) 



7*200 



while the characteristic density is given by 



201 



3 ln(l + c) - c/(l + c) ' 



(40) 



(41) 



Therefore, it remains to specify the concentration parame- 
ter as a function of the halo mass, which we have labeled 
M200 above to remind that it corresponds to the mass de- 
fined by the nonlinear density contrast S = 200. Although 
we shall also consider fits to numerical simulations that 
have already been proposed in the literature, we shall find 
in Sect. [5] below that a better match to numerical simu- 
lations for the density power spectrum (especially in the 
high-fc tail) is obtained with the following prescription, 



c(M 2 



10.04 



Mi 



2 x lO 12 ^- 1 ^ 



o 



(1 + *)- 



(42) 



which is in t ermed iate between the beha viors found in 
iDolag et all (|2004D and iDuffv et~aTI (|2008f ). However, this 
should not be trusted beyond the range studied in this pa- 
per (z < 3). 

Since we obtain Eq. (|42l) from a comparison with the 
density power spectrum measured in simulations, it can 
be understood as an "effective" concentration parameter 
that also includes the mean effect of halo substructures. 
Although it is possible to build more soph isticated halo 
mode ls that explicitly describe substructures (jGiocoli et al.l 
we do not investigate such refinements here. In the 
same spirit, it could be possible to include the mean effect 
of baryons on the total matter profile of the halos through 
such an "effective" concentration parameter (or another 
choice for the shape of the halo profile than the NFW for- 
mula ((391) ). However, this is also beyond the scope of the 
present study, and we leave such investigations to future 
works. In any case, it is clear that any modeling of halo 
profiles (cither at an "effective" level or through detailed 
explicit models) can be incorporated in the framework we 
describe in this paper. In particular, this could improve the 
accuracy at high k as simulations probe smaller scales and 
provide tighter (and more robust) constraints on halo prop- 
erties. 

For practical purposes it is convenient to use halo pro- 
files that have an explicit expression for the Fourier trans- 
form UM^k) defined in Eq. (|3~Tj) . especially if one is interested 
in high wavenumbers k where the integrand shows fast os- 
cillat ions. Let us recall that for the NFW profile (|3"9"|) one 
has (jScoccimarro et aDl200lh 



u M {k) = 



ln(l + c) - 



1 



— sin(cfcr s ) 
(1 + c)kr s 



■cos(kr s )[Ci[(l + c)kr s ]-Ci(kr s )] 
-sm(kr s )[Si[{l + c)kr s }-Si{kr s )]\, (43) 



where Ci(x) = — cost dt/t an d Si(a:) = smtdt/t 
are t he Cosine and Sine integrals (jAbramowitz k, Stegunl 
1970). 



4. N-body simulations 

In this section we describe the details of our N-body simula- 
tions, and we show the results of some tests of their possible 
systematic errors. 



4.1. Set up 

We adopt a flat ACDM 
mological parameters (57 m 
(0.279, 0.046035, 0.701, 0.817, 0.96), 
with 5 year observation of WMAP dKomatsu et al 



model with cos- 
Q h ,h,a s ,n s ) = 
which is consistent 



2009). 



We use a publicly available code, CAMB (jLewis et al. 2000), 
to compute the linear power spectrum. All the initial 
conditions for the simulations are created by adding 
displacements to particles put on a regul ar lattice by 
second-order Lagrangian perturbation theory (jScoccimarrol 
119981: ICrocce et ai1l2006l) at z = 99, and evolved b y a pub- 
licly available Tree Particle-Mesh code, GADGET2 (jSpringell 
l2005f ). We employ N = 2048 3 particles in four periodic 
boxes with different box sizes (I/box = 512, 1024, 2048 and 
4096/i _1 Mpc). We call these runs as L9-N11, L10-N11, 
Lll-Nll and L12-N11, respectively. The softening lengths 
for the tree force are set to be 5% of the mean inter-particle 
distance. 

In addition to the four main simulations, we run five 
smaller simulations (L9-N8, L9-N9, L9-N10, L10-N9 and 
L11-N10) for convergence tests. The box sizes and number 
of particles of these simulations are summarized in Tab. [TJ 



Table 1. List of N-body simulations presented in this pa- 
per. The box size, Lbox, is shown in units of /i _1 Mpc. 



name 




N L/iS 


Main 


Testl 


Test2 


L9-N11 


512 


2048 








L10-N11 


1024 


2048 


V 






Lll-Nll 


2048 


2048 








L12-N11 


4096 


2048 






V 


L9-N8 


512 


256 




V 




L9-N9 


512 


512 








L9-N10 


512 


1024 




V 




L10-N9 


1024 


512 








L11-N10 


2048 


1024 









4.2. Measuring power spectrum and two-point correlation 
function 

We basically follow the ordinary FFT method to measure 
both the power spectrum and the two-point correlation 
function. Before assigning particles to the mesh, however, 
we fold the particle distribution into a smaller box to avoid 
the systematic error caused by the assignment on a finite 
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Fig. 2. Wavenumber ranges covered by our N-body simulations, from the fundamental mode (left edge) to the mean 
inter- particle scale (solid right edge) and to the softening scale (dashed right edge). Left: the four main simulations, 
middle: simulations used in Test 1, right: simulations used in Test 2. 




1.04 



1.02 



3 1.00 

tL 



0.98 



0.96 



I I I I 
1 z=1 




1 L9-N10 


i _ 
i - 


L9-N9 


i - 
i - 


'_ L9-N8 


i 

/ _ 

1 — 


v ^\ 


i - 


s 

- \ 

\ 

- \ 

\ 


«. i ■— 

r-— ^ 

/ 

/ - 

/ 


- V 

i ,, i i •? 




0.1 0.2 0.5 1.0 2. 5. 10.0 

k [h Mpc" 1 ] 



0.1 0.2 0.5 1.0 2. 5. 10.0 0.1 0.2 0.5 1.0 2. 5. 10.0 

k [h Mpc" 1 ] k [h Mpc" 1 ] 



Fig. 3. The ratio of the power spectrum to the highest resolution run, L9-N11. The three lines correspond to L9-N10, 
L9-N9 and L9-N8 from top to bottom. Left: z = 3, middle: z = 1, right: z = 0.35. 



number of grid points, w hich is important on small scales 
near the grid interval (see lColombi et aTll2009L for a similar 
discussion) . Namely, we perform the following operation to 
the particle position, x: 

x^x%(L box /2"), (44) 

where the operation, a%b, gives the reminder of the di- 
vision of a by 6, and we vary the integer n from to 6. 
After that, we assign particles on the 2048 3 grid points 
of the folded small box with side length of Lbo x/2" using 
the CIC algorithm (jHocknev fc Eastwood! Il98lh . We then 
Fourier transform the density contrast and take the aver- 
age of |5(k)| 2 within fc-bins to obtain the power spectrum. 
The binning is chosen to be logarithmic on small scale with 
20 bins per decade (k > 0.3h Mpc" 1 ), but we adopt a lin- 
ear binning with the interval Ak = 0.005/iMpc" 1 on large 
scale (fc < 0.3h Mpc" 1 ) to see the feature of BAOs more 
clearly. While the folding procedure makes the systematic 
effect caused by assignment smaller, it reduces the number 
of available modes. We choose the value of n for the folding 
depending on the wavenumber so that the number of modes 



is large enough. We plot in Fig. [7] the resulting statistical 
error (1-er level) estimated from 



P(k) y/N modc (ky 

where N mol ±c(k) denotes the number of independent 
mode in that bin. T his is exact when <5(k) is Gaussian 
(jFeldman et al.l [1994), and this can be shown to re- 
main a good es t imate in the current situation (see e.g., 
Takahas hi et al.l (|2009() ). One may notice a dip at k ~ 
0.3/iMpc" 1 . This corresponds to the change of binning from 
linear to logarithmic. Another feature is a characteristic zig- 
zag pattern on smaller scale. This pattern corresponds to 
the change in the value of n used in the folding. Even after 
reducing the number of available modes by folding, how- 
ever, we keep enough Fourier modes to achieve accuracy of 
~ sub-percent level on small scale (k > 1/iMpc -1 ). 

The two-point correlation function is measured in a sim- 
ilar manner. We simply perform an inverse FFT to the three 
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Fig. 4. Same as Fig. [21 but for two-point correlation function. 
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Fig. 5. The ratio of the power spectrum to the largest volume run, L12-N11. The three symbols correspond to L11-N10 
(triangles), L10-N9 (circles) and L9-N8 (diamonds). Left: z = 3, middle: z = 1, right: z = 0.35. 



dimensional power spectrum calculated on grid points. We 
estimate the statistical error using the formula: 

|A£(z)| 2 = 64n 4 Jyk 2 {j (kx)P(k)} 2 , (46) 

where the quantity jo is the spherical Bessel function of 
the first kind. We use the power spectrum measured from 
the same simulatio n in the integrand for consistency, see 
iTaruva et all (|2009D for more details. 

4.3. Convergence tests 

One may naively think that one can model the power spec- 
trum by combining simulations with different box sizes. 
However, it is not trivial to determine which simulation 
gives the best measurement on a given scale. The simula- 
tion in a larger volume lacks the power on small scale, while 
that in a smaller box does so on large scale. 

A characteristic wavenumber for the finite box size is 
given by 2tt/L^ ox , the fundamental mode. For the finitc- 
ness of mass/force resolution, one may define two important 
wavenumbers. First, the initial density fluctuation cannot 
be input on scales smaller than the mean inter-particle dis- 
tance (27r/ibox x N 1 / 3 in wavenumber). Second, the force 
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Fig. 7. 1-er statistical error of the power spectrum (frac- 
tional) estimated from Eq. |g§) for L9-N11, L10-N11, Lll- 
Nll and L12-N11 (from top to bottom). 



calculation is not accurate below the softening scale. This 
wavenumber is given by 2n/L hox x iV 1 / 3 x 20 for our sim- 
ulations (i.e., 5% of the mean inter-particle distance). One 
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Fig. 6. Same as Fig. [5j but for the two-point correlation function. 



can basically trust the simulation results from the funda- 
mental mode to the inter-particle distance, or at most to 
the softening scale. These wavenumbers are shown in the 
left panel of Fig. [2] for the four main simulations. 

After nonlinear evolution, however, the initial small sys- 
tematic error on one scale may propagate to different scales 
due to the coupling between different modes. This makes 
the situation complicated. To separate the two systematic 
effects caused by the hniteness of mass/force resolution 
and simulated volume, we perform two tests using addi- 
tional simulations in this subsection. We will show that 
both hniteness effects lead to a lack of power on certain 
scales. We especially focus on the scales where these effects 
are more than 1%, and avoid to use the data points on these 
scales in later discussions. 

4.3.1. Test 1: finite mass/force resolution 

We first examine the effect of the hniteness of the 
mass/force resolution. We use L9-N8, L9-N9, L9-N10 and 
L9-N11 for this test (see Tab. [T|). These simulations have 
the same volume but different number of particles, so that 
we can see the systematic effect purely from the hniteness 
of the resolution (see the middle panel of Fig. We set 
the same random linear density held in creating the initial 
conditions of the four simulations to make the compari- 
son clearer. Since the three additional simulations, L9-N8, 
L9-N9 and L9-N10, have the same resolution as L12-N11, 
Lll-Nll and L10-N11, respectively, we can assess the im- 
pact of the systematic effect for the " main" simulations by 
this test. 

The results are shown in Fig. [3] We plot the ratio of 
the power spectrum measured from the three lower res- 
olution runs to that from the highest resolution one (i.e., 
L9-N11). One can clearly see that all curves blow up at large 
wavenumbers. This corresponds to the scale where the shot 
noise contribution becomes important. In addition to this 
feature, the curves show a lack of power on intermediate 
scales, which is more important for lower resolution runs at 
higher redshifts. We can interpret this as follows. At higher 
redshifts, the lack of power on around mean inter-particle 
scale survives. But the power generated by the pure non- 
linear growth gradually dominates over the initial power 



and the difference between the four runs become smaller at 
lower redshifts. 

We also test the effect of finite resolution on the two- 
point correlation function. Similar trends can be seen in 
Fig. H] for the two-point correlation function. 

4.3.2. Test 2: finite box size 

We next discuss the effect of the hniteness of the simulated 
volume. We use L9-N8, L10-N9, L11-N10 and L12-N11 for 
this test (see Tab. [I). These runs have the same mass/force 
resolution but different box size (see the right panel of 
Fig. [2]). Again, since the three smaller simulations have 
the same volume as the main runs, L9-N11, L10-N11 and 
Lll-Nll, we can test the possible systematics of these main 
simulations. 

We show the results in Fig. [5] In contrast to the test 
for the finite resolution, the data points seem noisy. This 
is because we cannot set the same random realization of 
the linear density field for simulations with different box 
sizes. When one sees the result of L9-N8, the systematic 
effect is prominent at 0.3 ~ 2/iMpc , and is growing with 
time. This is contrary to what we see in the finite resolution 
test, which as we showed was decaying. This growing nature 
of the systematics suggests that the effect comes from the 
mode coupling of large-scale modes. The systematic effect 
is not important (i.e., below 1% level) for the rest of the 
simulations. 

We test the convergence of the two-point correlation 
function which is shown in Fig. [5] Unfortunately, however, 
we cannot determine a strict range of the trustable scales 
due to large statistical error. Since the two-point corre- 
lation function is a weighted sum of the power spectrum 
over all scales, a large error bar at a certain wavenumber 
may propagate to the two-point correlation function over a 
wide range of scales. Furthermore, the two-point correlation 
function is not positive definite unlike the power spectrum. 
Thus, the size of its fractional error can be very large on 
scales where the two-point correlation function is close to 
zero. Finally, we cannot control the size of error bars by 
changing the width of the bins, since it does not explicitly 
depend on this width [see Eq. (|4"B])]. 

The ratio of the two-point correlation function on scales 
above > 5/i _1 Mpc seems consistent with unity at all the 
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Fig. 8. The probability F2n(l/k) that a Lagrangian pair of 
distance q = 1/k belongs to different halos, from Eq.([5T]). 
We show our results for z — 0.35, 1, and 3. 



three redshifts. As seen in Fig. [51 the finite resolution effect 
is important below < 5/i _1 Mpc for L9-N8. This prevents 
us from a clear test for the finite box size since the four sim- 
ulations used in this test have the same resolution as L9-N8. 
We conclude here that we do not find any clear evidence of 
systematic error on the two-point correlation function aris- 
ing from the finiteness of volume, and we will determine 
which simulation to use following only the resolution test. 
We leave more complete tests for the finite volume to future 
works. 



5. Comparison with numerical simulations 

We first show in Fig. [5] the quantity F 2 h(1/&), which is 
defined by Eq. (j2~lj) and enters the 2-halo term (|27l) , at red- 
shifts z = 0.35, 1, and 3. At fixed comoving scale, q — 1/k, 
it decreases at lower redshift, since the probability Fin for 
a Lagrangian pair of distance q to belong to the same col- 
lapsed halo grows with time, as larger scales turn nonlin- 
ear, and F2H = 1 — Fih- We can note that it decreases very 
slowly at high k. In fact, since we consider an initial power 
spectrum with a primordial slope n s = 0.96 that is smaller 
than unity, the rms linear density contrast cr(q) of Ea. dlOp 
does not diverge to infinity as q — >• but goes to a finite 
asymptote. This implies that the probability Fm{l/k) does 
not go to zero at high k (i.e. on small scales), but reaches the 

nonzero asymptote F2h(0) = ^ L ^" dvf[v)/v. Although 
this value is unlikely to be exact, because of the approxi- 
mations involved in the derivation of F2h(<z), it is possible 
that for such power spectra without significant power on 
small scales there always remains some fraction of mat- 
ter which has never experienced shell crossing and remains 
described by the single-stream fluid approximation. Even 
though this is an interesting point from a theoretical per- 
spective, regarding the dynamical properties of the system, 
this is beyond the subject of the present work and it plays 
no role for our purposes, since at high k the power spectrum 
is clearly dominated by the 1-halo term, as we shall check 
in Fig. below. On the other hand, we can see that if we 
require a high accuracy, the deviation of F2H from unity at 
k ~ 0.3/iMpc -1 for z = 1 cannot be completely neglected 
and has some effect on the 2-halo term (|27[) . 



We now compare our model for the density power spec- 
trum, that combines a systematic perturbative approach 
with a phenomenological halo model, with results obtained 
from numerical simulations. We show in Fig. [§] the power 
per logarithmic interval of fc, defined as 



A 2 (fc) = 47rfc 3 P(fc), 



(47) 



for the three redshifts z = 0.35, 1, and 3. The sudden rise of 
the N-body results at very high k (which is more apparent 
in the right panel at z = 3) is due to the finite resolution and 
shot noise, see Sect. 2) This marks the highest wavenumber 
where we can compare the simulations with our theoretical 
models. 

As discussed in Sect. 13.11 we can check that the 2-halo 
term A?, H remains well-behaved at high k and is subdom- 
inant with respect to the 1-halo term A 2 H . It falls some- 
what below the linear power A 2 , contrary to the one- loop 
resummation of A 2 ort given by the direct steepest-descent 
scheme , which becomes very close to A| as seen in lValageai 
(|2007aD , because of the prefactor F 2H (l/fc) in Eq. (f2T)) . 

In agreement with the discussion of Sect. 12.31 we can 
also check that the 1-halo term shows a fast decline at 
low k. At high k we can see that it is possible to reach 
a good agreement with the numerical simulations by using 
an appropriate prescription for the concentration parame- 
ter c(M), such as the one given in Eq. (|4"2")) and represented 
by the green solid line. The formula given in Eg. (1421) is 
obtained by looking for values of the free parameters (the 
normalization and the two exponents) that provide a rea- 
sonable match with the power spectrum measured in the 
simulations at high k (looking among a few values close to 
the fits already proposed in the literature). Nevertheless, it 
is interesting to note that using fits for c(M) proposed in 
previous analysis of the halo profiles formed in numerical 
simulations, one o btains predi ctions that are eit her larger 
(|Dolag et al.l l2004) or smaller (|Duffv et al.ll2008D than the 
one associated with Eq. (|4"2")) . This shows that the nonlin- 
ear power spectrum is fully consistent in this range with 
a simple halo model, such as (|35[) . and with the proper- 
ties of halos seen in numerical simulations. This also gives 
an estimate of the dependence of the power spectrum on 
the pr escription used for c(M) . In agreement with prev ious 
works (|Huffenberger fc SeliaklboollGiocoli et al.ll2010D . we 
recover the fact that the nonlinear power A 2 (fc) is larger 
and steeper at high k for concentration relations c(M ) that 
have a larger normalization and a steeper dependence on 
M. A nice feature is that the nonlinear power is largely in- 
dependent of the details of the halo model up to A 2 < 100, 
so that models such as the one studied in this article re- 
main quite predictive. Even at higher fc, we can see that up 
to A 2 < 10 3 , or fc < 100/iMpc -1 , using any of these pre- 
scriptions for c(M) provides a reasonable estimate of the 
power spectrum and even fares better than th e direct fit to 
P(k) that was obtained from older simulations (jSmith et al.l 
120031 ). This suggests that models based on phenomenolog- 
ical ingredients such as the halo model may prove more 
robust than direct fits to numerical results. However, our 
model for P(fc) should not be trusted beyond the domain 
where it has been checked, that is fc < 100/iMpc -1 and 
z < 3. 

In order to obtain predictions at much higher redshifts 
and wavenumbers, one should use a prescription for c{M) 
that is based on some physical arguments rather than sim- 
ple fits such as Eq. (j4*2"j) . This would probably lead to a loss 
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Fig. 9. The power per logarithmic interval of k, as defined in Eq. (|4"Tj) , at redshifts z — 0.35, 1, and 3. The symbols are 
the results from the numerical simulations described in Sect. HI The black solid line is the linear power spectrum, A? , 
and the blue dotted line that is somewhat below at high k is the 2-halo contribution (|27l) . A?, H . The steep solid line Af H 
is the 1-halo contribution (|3"5")) . for the halo model described in Sect. 13.21 The green solid line is the full nonlinear power 
spectrum, A 2 = A| H + A 2 H . The dashed (resp. dotted) green line, that is slightly larger ( resp. s maller) at high k, is 
the result obtained by using the concentration relation for c(M) given bv lDolag et ail ([2004) (resp. iDuffv et al.1 (l2008ln 
instead of Eq. (|4"2")l . The mag enta dot-dashed lin e A|, that is somewhat below these new N-body results at high fc, is the 
fit to older simulations from lSmith et all (|2003[ ). 



of accuracy in the range tested in Fig. [9j as compared with 
the use of Eg. (14"2"1) . but this is likely to be more robust as 
we extrapolate to other regimes. However, we shall not in- 
vestigate the building of physical models for c(M) in this 
article, as this is a topic by itself. 

We can see in Fig. [HI that the full nonlinear power spec- 
trum (|22[) obtained by our approach, combining a pertur- 
bative expansion with a phenomenological halo model, is 
able to reproduce the results measured in numerical simu- 
lations, up to a reasonable accuracy. Let us point out that 
the perturbative part, P pe rt(k), which dominates on large 
scales, contains no free parameter, since it is given by a sys- 
tematic perturbation theory. The 1-halo contribution that 
dominates on small scales clearly contains some parame- 
ters, through the choice of the halo mass function and halo 
density profile. However, the mass function and the shape 
of the halo profile are already set by other measures from 
numerical simulations, so that the main free parameter is 
the concentration c(M), which is not as well constrained. 
However, as seen in Fig.[3J this uncertainty only affects the 
very high-fc tail. Moreover, even in this region the resulting 
model is competitive with direct fits to the power spectrum 
measured in older simulations. Therefore, models such as 
the ones developed in this article should prove useful to 
obtain reliable predictions for P(k). Another advantage of 
such models, as compared with simpler fitting formulas, is 
that they should be quite accurate on large scales since they 
are consistent with perturbation theory, up to the order of 
truncation of the computation. 

Thus, we plot in Fig. [To] the ratio of the power spec- 
trum to a "no- wiggle" linear power spectrum Pz, s (fc) which 
does not include baryonic acoustic oscillations, in order to 
see more clearly the low-fc behavior. We can check that in 
this range our model, based on the direct steepest-descent 
resummation, shows a good match with the numerical sim- 
ulations. As could be exp e cted, it fares better than a di- 
rect fit from ISmith et al.l (|2003|) which was not designed 



to model this range with a high accuracy (but still re- 
mains surprisingly good). Furthermore, the use of the di- 
rect steepest-descent resummation proves to provide sig- 
nificantly more accurate results than the standard pertur- 
bation theory, truncated at one-loop order, given by the 
well-known expressions 

Aioop (k) = Pl (A:) + P22 (k) + P31 (k) , (48) 

where formula s for the terms P22 and P31 may be found 
for instance in iBernardeau et al.l (|2002f ). Let us recall here 
that the term P pert (fc) given by this direct steepest-descent 
scheme contains no free parameter, nor any interpolation 
procedure, and is consistent with standard perturbation 
theory up to one-loop order (i.e. the difference between 
-Piioop and Pp Cr t is due to the partial resummation of higher 
order terms). We can see in the left panel, at z = 0.35, that 
around k ~ 0.14/iMpc" 1 the curve P per t is slightly above 
the full nonlinear power spectrum P. This is due to the pref- 
actor P2h(1/A;) in Eq. (f2~Tj) , which is slightly below unity. 
However, this is only a very small effect on these scales. 
On the other hand, at higher k (e.g., k > 0.26/iMpc _1 at 
z = 0.35) the nonlinear power spectrum rises above P pe rt 
and keeps growing, while P per t remains close to Pl at high 
k as seen in Fig. [9l This is due to the 1-halo contribu- 
tion, which starts being non- negligible. However, on these 
scales the dependence on the details of the halo model is 
extremely weak. Indeed, the three green curves plotted in 
Fig. [SI associated with the prescription (l4"2l for the c oncen- 
trati on c(M) and the two fits given by iDolag et al.l (|2004) 
and lDuffv et al.l (|2008l) . are also plotted in Fig.[10l However, 
they almost exactly fall on top of each other and cannot be 
distinguished in the figure. 

Thus, Figs. l9l and [TOl show that by combining perturba- 
tion theories and halo models it is possible to obtain a good 
model for the nonlinear density power spectrum, both on 
quasi-linear and highly nonlinear scales. However, we can 
see in Fig.|5]that in the intermediate regime, where A 2 ~ 5, 
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Fig. 10. The ratio of the nonlinear power spectr um P(k) to a smooth linear power spectrum Pl s without acoustic 
baryonic oscillations, from lEisenstein fc Hul (|l999h . The points with error bars are the results from N-body simulations. 
The black solid line Pl is the linear power spectrum, and the upper blue dashed line Piioop is the prediction of standard 
perturbation theory up to one- loop order. The green solid line P is the full nonlinear power spectrum, from Eq. (|22[) . 
wher eas the lower blue d otted line P per t is the perturbative part. The magenta dot-dashed line P$ is the fit to simulations 
from lSmith et al.l (|2003l ) (which was not specifically designed to reproduce the baryon oscillations). 



our predictions fall below the N-body results. This is also 
apparent in the high-fc parts of Fig. 1 101 where the full non- 
linear prediction P(k) starts to grow more slowly than the 
power measured in the numerical simulations. This regime 
corresponds to the transition between the 2-halo and 1-halo 
contributions (see Fig. [9]) and as such it is at the limit of 
validity of the approximations used for both terms. 

On the perturbative side, that is the 2-halo term, the 
discrepancy can be due to the truncation at one-loop order 
of the perturbative term. Indeed, we can expect that by 
going to higher orders we can extend the range of validity 
of the the perturbative term Pp Cr t and push the downturn 
shown by the blue dotted curve in Fig. [10] to higher k. 
Within such resummation schemes this means that we in- 
clude all diagrams up to n loops, and partial resummations 
for higher order terms. We can note that the discrepancy 
looks somewhat more severe at z = 3 in Fig. ® More pre- 
cisely, the range of k where there is a noticeable mismatch 
before the 1-halo term becomes dominant (larger than Pl) 
is somewhat more extended than at z = 0.35. This agrees 
with the results of IValageasl (|2010al ). where it was found 
(within the Zeldovich framework) that the scope of pertur- 
bation theory is somewhat greater at higher redshift for 
CDM power spectra, in the sense that the range where 
higher order perturbative terms are important (i.e. larger 
than the non-perturbative correction associated with shell 
crossing effects) is wider and that the perturbative expan- 
sion makes sense up to higher orders. However, we shall not 
try to go beyond one-loop order in this paper and we leave 
such a task to future works. 

On the non-perturbative side, it is clear that by defi- 
nition halo models are only phenomenological models and 
do not provide a systematic tool. Even if we admit that it 
makes sense to "recognize" halos in the density field (even 
though this is a descriptive tool rather than a direct ingre- 
dient of the equations of motion), it is clear that realistic 
halos are not spherical nor fully virialized. In particular, 
apart from the inaccuracies introduced by the spherical ap- 
proximations, outer halo radii are not fully virialized and 
particles do not instantly loose all memory of their initial 
conditions, so that the average in Ea. (|33[) is only approxi- 
mate. As we have discussed in Sect. 12.31 such approxima- 



tions also prevent us from exactly satisfying momentum 
conservation. Since these effects mostly apply to the outer 
parts of the halos they can be expected to have an impact 
on the moderately low-fc part of the power spectrum, that 
is the transition region. Another effect that is likely to play 
a role is the truncation of the halo profiles at the radius 
r2oo- This may be investigated in a rather straightforward 
manner by studying the change in P(k) as we modify this 
density threshold, and we shall consider in Sect. 16. ll below 
the results obtained for 5 — 50. 

On a more fundamental level, since the splitting be- 
tween 2-halo and 1-halo contributions is necessarily approx- 
imate (it cannot be derived in a systematic fashion from 
the equations of motion and always involves some approx- 
imations, at least within analytical frameworks) it is not 
surprising that the match is not perfect, especially in the 
transition range. 

Another approximation used in this work is the assump- 
tion of smooth halos, defined by the regular profile (131?)) . As 
noticed in Sect. 13.21 it is pos sible to extend the h alo model 
by including substructures ()Giocoli et al.l I2010T ) , however 
this only has an effect on the very high-fc tail of the power 
spectrum (as could be expected since these are small-scale 
modifications). We do not investigate such modifications 
in this work, as the agreement we obtain with numerical 
simulations is already reasonably good on these scales and 
we prefer not to introduce additional parameters. These 
effects are also degenerate with the mean concentration re- 
lation c(M), as seen in Fig. Oby the comparison between 
the dashed and dotted green curves , assoc iated with differ- 
ent pr escriptions from iDolag et al.l (|2004l ) and lDuffv et al.l 
(2008J). In particular, as pointed out in Sect. 13.21 the pre- 
scription (I42[) that we have obtained by requiring a good 
match at high k can also be seen as an effective model, that 
describes the combined effects of both the concentration of 
the mean radial profile and the presence of small substruc- 
tures. At some stage, if one wants to build a model where 
the halo properties are independently obtained from the 
measure of individual halos in N-body simulations, so that 
c(M) can no longer be tuned, it may nevertheless prove 
necessary to explicitly include such effects. 
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6. Dependence on various ingredients of the model 

We now investigate the impact of various ingredients of the 
model on the predictions obtained for the density power 
spectrum. 

6.1. Halos defined by 5 = 50 

As we have discussed in the previous section, the lack of 
power on intermediate scales seen in Fig. [SJ around the 
transition between the 2-halo and 1-halo terms, may be due 
in part to the truncation of the halo profiles at the radius 
^200- Indeed, realistic halos extend to larger radii, which 
should generate some extra power on large scales. In order 
to take into account such outer radii, while keeping the 
mass M that enters the definition of the halo mass function 
as the halo mass within the truncation radius (so as to 
avoid overcounting and violations of mass normalizations) , 
we consider defining halos by a smaller nonlinear density 
contrast, such as 6 = 50. This lower threshold means that 
we include outer shells (as compared with the case 5 = 200), 
while halos still remain well separated from the background. 

For the halo mass function we still use the scaling func- 
tion (|38l) . but the linear density contrast is now 6l = 
.F _1 (50), which gives for instance 6l — 1-50 at z = 1. This 
automatically satisfies the normalization (IT51) a s well a s the 
large- mass tail (fT^)l . It was shown in I Valageasl (|2009al) that 
this recipe provides a reasonable match to numerical sim- 
ulations (but not as good as for S = 200) for halos defined 
by a density threshold 5 = 100, and we shall assume that 
it still provides a reasonable model for S = 50. This should 
be sufficient for our purposes, as the 1-halo term needed 
for the density power spectrum only depends on integrals 
over the halo mass function and we have seen that both the 
normalization and the large-mass tail are correctly repro- 
duced. 

We still use the NFW profile of Ea.p?]). but we now 
need to obtain the scale radius r s in terms of r 50 . To do so, 
we now define the concentration parameter as 

C (M 50 ) = (49) 
r s 

and the halo characteristic density now reads as 

- 51 c 3 

3 ln(l + c) — c/(l + c) 

We again look for a prescription for c(M$q) of the form of 
Eq. (|4"2"]) . and we find in Fig. [TT] below that a reasonable 
agreement with the numerical simulations can be achieved 
by using 

This is mostly set by the behavior of the power spectrum 
at high k, where A 2 > 100, as seen in Fig. [5] from the 
comparison between various models for c(M) for halos de- 
fined by S = 200. Thus, the properties of the density power 
spectrum on large scales, where A 2 < 10, are almost inde- 
pendent of the choice of parameters in Eq. ([51]) (restricted 
to a reasonable range). 

We compare in Fig. [11] the power spectra obtained on 
nonlinear scales by using halos defined either by 6 — 200 or 




k [h Mpc" 1 ] 

Fig. 11. The power per logarithmic interval of fc, as in 
Fig. [§J but for halos that are defined either by a density 
contrast S = 200 (solid curves, identical to those of Fig. [3]), 
or by 5 — 50 (dot-dashed curves). Both the 1-halo contri- 
bution Am and the full nonlinear power spectrum A 2 are 
shown. 
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Fig. 12. The ratio of the nonlinear power spectrum P(k) 
to a smooth linear power spectrum P^ s without acoustic 
baryonic oscillations, as in Fig. [12J for halos defined by 
5 = 200 (green solid curves, identical to those of Fig. [TU|) . 
or by 6 = 50 (red dot-dashed curves). 



S = 50 (while the perturbative term P per t(fc) is not modi- 
fied). As expected, truncating halo profiles at a larger radius 
yields some extra power at large k and improves somewhat 
the agreement with numerical simulations in the intermedi- 
ate range, where A 2 ~ 5. However, the change is quite small 
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and the match to the N-body results is not perfect yet. In 
fact, at redshift z = 0.35, we obtain a smaller power spec- 
trum around A 2 ~ 100 by using halos defined by 5 = 50. 
However, one may modify this by using a different concen- 
tration relation c(Mso) at that redshift. In any case, these 
results suggest that the truncation at r2oo is not the unique 
reason for the discrepancy between the model and N-body 
simulations in this transitory range, and as explained in 
Sect. [5] another possible source of inaccuracy is the partial 
resummation of high-order terms in the perturbative term 

P pert • 

Next, we compare in Fig.[T2]the power spectra obtained 
on quasi-linear scales. We can see that on these very large 
scales, k < 0.3/iMpc -1 , the modification of the 1-halo term 
has a very weak effect. However, it improves somewhat the 
agreement with the simulations at low redshift, z = 1 (at 
k ~ 0.3/iMpc -1 , as the gap below the N-body results in the 
intermediate range is repelled to slightly higher k), although 
it leads to a slight overestimate for P{k) at z = 0.35. On 
the other hand, it appears to make almost no change over 
these scales at higher redshift, z = 3. In ag r eement with 
the discussion in Sect. [5] and with IValageaa (|2010alh this 
expresses the fact that for CDM power spectra the scope of 
perturbative expansions is somewhat broader at higher z, in 
the sense that the range of k where higher order terms play 
a role is wider. This corresponds to the interval [£;n oop , fc s . c .] 
where terms beyond linear order are already significant, as 
compared with the linear power Pl, but still larger than 
the non-perturbative correction associated for instance with 
shell crossing effects. 

That our results for the density power spectrum only 
show a very weak dependence on the details of the halo 
definition (by 6 — 200 or 50) can actually be considered 
as a reassuring sign, rather than a problem in the attempt 
to improve the agreement with simulations. Indeed, if the 
splitting (f2"2"j) and all subsequent computations were exact, 
halos would only appear as an intermediate tool and would 
disappear in the final result, since the matter density field 
is independent of how we split it over halos at later stages. 
Therefore, the final result for P{k) should not depend on 
the details of the definition of the halos. In practice, this 
cannot be the case because the splitting (|2"21 itself is only 
approximate as long as we consider spherical halos, and 
the subsequent treatment of both 2-halo and 1-halo contri- 
butions involves some approximations, as in Eqs. (|27p and 
(|3"2"|) . This leads to some artificial dependence on the de- 
tails of the halo model, but as shown in Figs. [TT] and [12] 
this effect is rather small and argues for the robustness of 
the model. 

Therefore, Figs. [TT] and [T3] suggest that there is still 
room for systematic improvement, by including higher or- 
der terms in the perturbative contribution P pcr t. However, 
this may come at the price of heavier and slower compu- 
tations, which we do not investigate here. On the other 
hand, they show that predictions on large scales are largely 
independent of the details of the underlying halo model. 

6.2. Impact of the "1-halo" counterterm 

We have seen in Sects . [51 and 16.11 that quantitative details of 
the 1-halo term can have some effect on the high-fc tail (the 
concentration c(M) of the halo profiles) or on the interme- 
diate range where A 2 (fc) ~ 5 (the truncation radius of the 
halos). We investigate here the role of a more fundamental 
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Fig. 13. The ratio of the nonlinear power spectrum P(k) to 
a "no- wiggle" linear power spectrum Pj, s (fc), as in Fig. [TD1 
at redshifts z — 0.35, 1, and 3. The green solid lines corre- 
spond to our fiducial model, with halos defined by a density 
contrast S — 200, and are identical to the green solid lines 
shown in Fig. [JTJ while the red dashed lines are obtained by 
setting to zero the 1-halo counterterm W(kqM) 2 of Ea. ([35|) . 



ingredient of the 1-halo contribution, namely the countert- 
erm W(kqM) 2 in Ea. ([35l) . Thus, focusing on the case of 
halos defined by the nonlinear density contrast <5 = 200, 
we compare in Fig. [TS] the predictions of our fiducial model 
with those that are obtained by setting to zero the 1-halo 
counterterm W(kqM) 2 of Eg. (1331) . Thus, we can see that, 
even though the 1-halo counterterm plays a negligible role 
at high k, it cannot be discarded at low k if we require a high 
accuracy. As explained in Sect. 12.31 this is also related to the 
fact that this counterterm is required to obtain the k 2 tail 
at low k expected for the 1-halo contribution (if we do not 
insist on momentum conservation). Discarding this coun- 
terterm leads to a nonzero asymptote at low k for Pm(fc), 
which eventually dominates over the linear power spectrum, 
which roughly decreases as Pl(A;) ~ k for CDM cosmolo- 
gies. On the scales probed in Fig. [TBJ k ~ 0.2/iMpc" 1 , we 
have not entered yet the very low k regime where this would 
give rise to a spurious divergence. However, we can clearly 
see the spurious extra power that would be generated by 
the lack of this counterterm. Note that on these scales there 
are no fitting parameters that could be tuned to compen- 
sate for this extra power. Indeed, the details of the halo 
model, such as the concentration c(M) (see Fig. [TO"]) , or 
the truncation radius of the profiles (see Fig. [T2"]) . only 
have a weak impact in this domain and could not bal- 
ance this extra power. On the other hand, the 2-halo term 
is also highly constrained, as the blue dashed and dotted 
curves in Fig.QJJ associated with standard 1-loop perturba- 
tion theory and steepest-descent resummation, give an es- 
timate of the possible scatter between the predictions that 
can be obtained from various perturbative schemes. This 
is most clearly seen around the peak at k ~ 0.13/iMpc -1 
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for z = 0.35 (as could be expected the extra power asso- 
ciated with such an incorrect modelization of the 1-halo 
term grows at lower z as larger halos form). By contrast, 
the reasonably good agreement with simulations obtained 
by our complete model (solid lines) suggests that it does 
not miss important physical ingredients (as discussed in 
Sect. [5] difficulties mostly arise at higher k, associated with 
the transition between the 2-halo and 1-halo contributions). 

6.3. Benefit of higher order perturbative terms 
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Fig. 14. The ratio of the nonlinear power spectrum P(k) to 
a "no- wiggle" linear power spectrum Pi s (fc), as in Fig. 1101 
at redshifts z = 0.35, 1, and 3. The green solid lines "resum" 
correspond to our model as in Fig. [TU1 whereas the blue 
dashed lines correspond to the approximation P2H = Pl 
(it cannot be distinguished from the linear power Pl on 
these scales at z = 3). 



Before we investigate other choices for the perturbative 
term P por t in the next sections, we first compare in Figs. [T4l 
and !15l our model with the more standard approach, where 
one simply uses the linear power Pl for the 2-halo con- 
tribution. Thus, we replace Eq. (|27l) by the approximation 
P-m{k) = Pi(fc), while we keep the same 1-halo contri- 
bution (|35[) . with its new counterterm that ensures a sat- 
isfactory behavior at low k. In agreement with previous 
works, we can see in Fig. Q3] that on the large scales as- 
sociated with the baryon acoustic oscillations the approx- 
imation P2B.{k) — Pt,(fc) is not sufficient to obtain a good 
match to the numerical simulations, contrary to the use 
of the one- loop perturbative term P pelt (k) obtained by the 
"steepest-descent" scheme recalled in Sect. 13. ll This means 
that on these scales the departure from the linear power 
is not yet due to the non-perturbative contribution associ- 
ated with the 1-halo term, but to higher order perturbativ e 
terms. This agrees with the results of Valageas] (|2010alh 
based on the Zeldovich dynamics, which show that many 
orders of perturbation theory are relevant before the power 
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Fig. 15. The power per logarithmic interval of k, as in 
Fig. [HI at redshifts z = 0.35, 1, and 3. We plot the curves 
obtained with our model (green solid lines) and with the 
approximation P 2 h = Pl (blue dashed lines), as in Fig.[T4l 



spectrum is dominated by the non-perturbative contribu- 
tion associated with shell crossing or halo formation (typi- 
cally one can go up to order P® at z = and Pf 6 at z = 3 
for a ACDM cosmology, at least for this simpler case). This 
motivates the use of perturbative approaches that include 
higher order contributions. 

On smaller scales, shown in Fig. [151 the difference be- 
tween our model for P2h(&) and the popular approxima- 
tion P2n(k) = PL(k) is negligible. This agrees with Fig. 
where we could see that at high k the one-loop resummed 
contribution to P2n(k) is near the linear power Pj,(fc) and 
subdominant as compared with the 1-halo contribution. 

6.4. Eulerian perturbation theories 

We consider here two other choices for the perturba- 
tive term P per t of the 2-halo contribution ([2~T]) . A first 
choice is simply to use the standard perturbation the- 
ory and to write at one-loop order P pcr t = Puoop(fc)) 
where Pi\ oop {k) was given in Eq. (j48|) . A second choice is 
to keep the diagrammatic expression given by the first line 
in Fig. [TJ but to replace the external double lines by the 
Gaussian-decay respon s e functi on (propagator) obtained by 
ICrocce fc Scoccimarrol (|2006bl lah instead of using the re- 
sponse function obtained by the steepest descent resum- 
mation scheme. This is not identical to t he "R PT" resum- 
mation used in ICrocce fc Scoccimarrol (|2008lh where the 
"bubble" in the upper right diagram of Fig. [T] is also re- 
placed by resummed two-point correlations, obtained from 
this same nonlinear propagator (which corresponds to the 
partial resummation of further diagrams). Here we keep 
linear two-point correlation functions in this "bubble" to 
keep factorizable integrals, as in Eq. (IA.29[) . since our goal 
is simply to estimate the dependence on the choice of re- 
summation scheme. Moreover, it is interesting to evaluate 
various factorizable schemes of this kind as they allow faster 



P. Valageas and T. Nishimichi: Combining perturbation theories with halo models 



19 




0.9 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

0.1 0.2 0.3 

k [h Mpc- 1 ] 



Fig. 16. The ratio of the nonlinear power spectrum P{k) 
to a "no-wiggle" linear power spectrum Pj, s (fc), as in 
Fig. [TU1 at redshifts z — 0.35, 1, and 3. The red dashed 
lines "SPT" correspond to the use of 1-loop "Standard 
Perturbation Theory" for the term P per t in the 2-halo con- 
tribution Q27p. whereas the blue solid lines "DR" corre- 
spond to the use of the "Decay ing Response" function from 
ICrocce fe Scoccimarrol (|2006aD in the external double lines 
of the two upper diagrams of Fig. [TJ The 1-halo term is the 
same as the one used in Sect. [5] 

comput ations. Two other differences with the prescription 
used in ICrocce fc Scoccimarrol (|2008t ) are the presence of 
the prefactor F<m{\/k) in Eq. (f2"T|) , and the fact that we 
do not introduce an extra tuning in the decay of the re- 
sponse function by replacing the linear velocity dispersion 
a v , which governs its high- A; cutoff, by a nonlinear veloc- 
ity dispersion obtained from fits to the nonlinear power 
spectrum from numerical simulations. Indeed, one of the 
main points of our approach is precisely that the term 
-Ppert (fc) in the 2-halo contribution (|27p should be obtained 
from systematic perturbation theories, to ensure a good 
accuracy and robustness for a variety of cosmologies, so 
that we avoid introducing fitting parameters in this factor. 
Then, these two alternatives, "standard perturbation the- 
ory" (SPT) and "Gaussian-like decaying response function" 
(DR), and the complete steepest-descent resummation used 
in the previous sections, agree with each other up to one- 
loop order. 

We first show in Fig. [TBI the power spectra obtained on 
large scales for redshifts z — 0.35, 1, and 3. We plot the re- 
sults obtained by using either the 1-loop "standard pertur- 
bation theory" , or the "decay ing response" function from 
ICrocce fc Scoccimarrol (|2006al ) in the external double lines 
of the two upper diagrams of Fig. [TJ for the term P per t 
in the 2-halo contribution (|2T|) . In both cases the 1-halo 
term is the same as the one used in Sect. [5j so that these 
curves only differ from those shown in Fig. [TJJ through the 
term P per t- As was already n oticed in Fig. [TO! and in agree- 
ment with previous works (ICrocce fc Scoccimarrol 120081: 
iTaruva fc Hiramatsdl20M ICarlson et al.ll2009t) . the 1-loop 



standard perturbation theory overestimates the power on 
these scales, especially at low redshift. The dashed curves 
shown in Fig. [16] are not identical to those shown in 
Fig. 1101 and in other works, because the perturbative term 
Ppert is multiplied by the prefactor F 2 n(l/k) in Eq. (|2"T|) . 
Moreover, the results also include the (small) 1-halo con- 
tribution. However, on these scales the prefactor P2h(1/^) 
remains very close to unity (see Fig. [5]) and it is not suffi- 
cient to significantly improve the agreement with the N- 
body results. The use of the "decaying response" func- 
tion leads to results that are very close to those obtained 
with standard perturbation theory and also overestimates 
the power spectrum on these large scales. Using the full 
"renormalized perturbation theory" (RPT) described in 
ICrocce fe Scoccimarrol (|2008h leads to a smaller prediction 
for P(fc). This can be understood from the fact that in- 
serting the "decaying response" function into the bubble of 
the upper right diagram of Fig. [TJ as in "RPT" , clearly sup- 
presses the nonlinear power and explicit computations show 
that this yields a good match to results from sim ulations 
(|Crocce fe Scoccimarroll2008l : ICarlson et al.ll2009f ). 
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Fig. 17. The power per logarithmic interval of fc, as in 
Fig. [9j at redshifts z = 0.35, 1, and 3. We plot the curves 
obtained using 1-loop "standard perturbation theory" (red 
dashed lines for the full nonlinear power A 2 , and lower red 
dotted lines for the associated 2-ha lo contribution) or the 
"decay ing response" function from ICrocce fc Scoccimarrol 
(|2006al ) (blue solid lines for the full nonlinear power A 2 , 
and lower blue dotted lines for the associated 2-halo con- 
tribution) . Both the full and 2-halo power spectra obtained 
with "SPT" are larger than their "DR" counterparts. 

Next, we show in Fig. [T7] the power per logarithmic 
interval of fc, as in Fig. [S] We plot both the full nonlin- 
ear power A 2 , as in Fig. [TBI an d the 2-halo contribution, 
for the "SPT" and "DR" cases. As already pointed out in 
Sect. 13.11 we can see that standard perturbation theory fails 
at a qualitative level as the 2-halo contribution keeps grow- 
ing at high fc and becomes very large, whereas the term 
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Ppert being associated with particle pairs that have not col- 
lapsed within a single object it should remain near Pl at 
most. Moreover, this spurious growth is not suppressed by 
the prefactor P 2 H(l/fc), which does not decrease sufficiently 
fast at high fc, as seen in Fig. [31 Even though this 2-halo 
contribution remains smaller than the 1-halo contribution 
on these scales, it leads to a significant extra power that 
worsens the agreement with the numerical simulations. As 
can be seen from Figs. l9l andfTTl this extra power cannot 
be compensated by modifications to the underlying halo 
model (such as the concentration c(M) or the halo trunca- 
tion radius), especially in the range 5 < A 2 < 200 where the 
dependence on the details of the halo model is rather weak. 
Of course, this problem would become increasingly severe as 
one include s higher order terms that grow increasingly fast 
at hig h k (|Bernardeau et al] 120021 : ICrocce fc Scoccimarrol 
2006)3). This means that standard perturbation theory can- 
not be used to build systematic models for the density 
power spectrum, such as those studied in this paper. As 
pointed out in Sect. 13.11 this provides a further motivation 
to investigate resummation schemes that would be better 
behaved. 

We can see that using the "decaying response" func- 
tion in the external double lines of Fig. [1] leads to a 2- 
halo contribution that remains small and close to the lin- 
ear power spectrum at high fc, as for the direct steepest- 
descent scheme used in Fig. |H1 This shows again the benefit 
of using resummation schemes instead of standard pertur- 
bation theory to obtain a well-behaved perturbative term 
-Ppert (k) that also agrees with standard perturbation the- 
ory up to the order of truncation. In fact, it happens that 
using this "decaying response" function yields a somewhat 
larger perturbative term than the one obtained in Sect. [5] 
on large scales, as seen in Fig. \TE\ and this extra power 
actually improves the agreement with the N-body results 
in the range where 1 < A 2 < 100. Indeed, as discussed in 
Sect.0 the implementation used in Fig. [5] underestimates 
the nonlinear power spectrum in this range, and the ad- 
ditional power due to the use of this "decaying response" 
function helps to bridge the gap. It actually provides a very 
good quantitative match, although at z = 3 we again un- 
derestimate somewhat the power spectrum in the transition 
range. Thus, for practical quantitative purposes this may 
provide a good variant of the model developed in this pa- 
per to compute the power spectrum, if one is interested 
in these transition scales. However, because on large scales 
this prescription actually gives a worse match to numerical 
simulations, as seen in Fig. 1161 the better agreement on the 
transition scales is probably not more than a lucky coinci- 
dence. Thus, it is likely to compensate a meaningful lack 
of power (due to the various reasons discussed in Sect. [5j 
such as the neglect of some higher order perturbative terms 
or the limitations of the description in terms of spherical 
virialized halos) by an unrelated artificial overestimate. 

6.5. Lagrangian perturbation theories 

The Lagrangian re-interpretation of the halo model devel- 
oped in this paper, leading to the expression (|2"o]l for the 
2-halo contribution, would most naturally lead us to con- 
sider Lagrangian perturbation theories. It is possible to 
substitute such schemes into Eq. (|2"6")) , which allows us to 
integrate in a self-consistent fashion over the q— dependent 
prefactor F2n(q). Unfortunately, this does not yield very 



successful results while leading to somewhat more complex 
expressions than those encountered in Eulerian frameworks. 
Therefore, in this section we illustrate the problems faced 
by current Lagrangian schemes by using the simpler expres- 
sion (|2"T[) , where the term P 2 h (<?) has been factorized as the 
prefactor P 2 H(l/fc)- Thus, as for the Eulerian perturbation 
theories used in the previous sections, we only need the pre- 
diction for the resummed power P pe rt (fc) to compute the 2- 
halo contribution. This simplifies matters and only makes 
a negligible change on large scales where P 2 h — 1. 
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Fig. 18. The ratio of the nonlinear power spectrum P(k) to 
a "no- wiggle" linear power spectrum Pj, s (fc), as in Fig. [TU1 
at redshifts z = 0.35, 1, and 3. The red dot-dashed lines 
"-FWt" correspond to the use of the Lagrangian-based 1- 
loop expansion (1531) for the term P per t in the 2-halo contri- 
bution (|2"7|) . whereas the blue solid lines "P m z" correspond 
to the use of the "modified Zeldovich" expansion (j56|) . The 
1-halo term is the same as the one used in Sect. [5j 

We consider two schemes that are inspired by a 
Lagrangian framework . A fir st approach is the one pre- 
sented in iMatsubara (|2008F) . which for our purposes 
amounts to expand P(k) after factorization of a Gaussian 

dampin g term e~ k (see also ICrocce fc Scoccimarrol 
(|2006bl ): IValaeeasI (|2010al) V where o v is the rms 
one-dimensional linear displacement (or the rms one- 
dimensional linear velocity, up to a time-dependent factor) 
given by 



47T 

T 



dfcP L (fc). 



(52) 



Thus, at 1-loop o rder the density power spectrum reads as 
(|Matsubarall2008r i 



P M at(fc) = e" fc ^ (P L + P 22 + P 31 + fc 2 a 2 Pi) , 



(53) 



where P 22 and P31 are the 1-loop terms obtained in 
the standard Eulerian perturbation theory, as in Eq.(|48p. 
Of course, Ea. (|53| could be obtained at once by requir- 
ing consistency with Eq. (l4*5|) at order P|. As discussed 
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Fig. 19. The real-space two-point correlation function at redshifts z — 0.35, 1, and 3. The symbols are the results 
from the numerical simulations described in Sect. 01 The black dashed line is the linear correlation, and the blue 
dotted line which is somewhat below at low x is the 2-halo contribution ^2H- The steep solid line xira is the 1-halo 
contribution, for the halo model described in Sect. 13.21 The green solid line is the full nonlinear correlation function, 

£, = + £lH- 



in IValaeeasI (|2007bl . 120081 l2010al) . and as can be seen 
from the computation of Lagrangian response functions 
(jBernardeau fc Valageas1l2008ll2010al ). the prefactor e'^^ 
is associated with the coherent displacement of density 
structures by the long wavelengths of the velocity field and 
should only appear in different-time quantities, under the 
form e~( Dl ~ D2 > k CT «o/ 2 ; where Di is the linear growth fac- 
tor at redshift Zi and a v o is the linear dispersion today, 
where D(z = 0) = 1. Thus, looking for an expansion un- 
der the form (|53p leads to an artificially strong Gaussian 
damping at high k, which has no physical meaning. In or- 
der to remedy to this problem, which can also be seen as a 
consequence of the truncation to a finite number of terms 
between the parenthesis in Eq. (|53[) (once we insist on look- 
ing for an expansion of this form) , we note that the nonlin- 
ear power spect rum generated by the Zeldovich dynamics 
(jZeldovichlll97(ih can also be written under this form, 

oo 

P z (fc) = e-^£pW(fc), (54) 

71=1 

where each term P^'Qc) scales as (Pl)™- Moreover, both 
the full nonlinear expression of the Zeldovich power spec- 
trum, Pz(k), (i.e. the resummation of the series ([5"4"))). and 
the expressi on of the terms (Pl)" of all orders, are explic- 
itly kn own (|Crocce fc Scoccimarrdl2"o""06bt IValageas 2007b, 
l2010al) . Then, a possible recipe is to use for the nonlin- 
ear power spectrum (more precisely the perturbative term 
Ppert) the "modified Zeldovich" expression 

Pmz(fc) = e- k2 °* ^(fe)+P (2) (fc) + £pi n) (fc)J (55) 

= P z (A) + e~ k ^ (P< 2 > (k) - P< 2) (*)) , (56) 

where P^ n > (k) is the term of order (Pl)" in the expansion of 
the form (|54[) for the gravitational power spectrum, whence 
P(2)(fc) = P 22 (fc) + P 31 (fc) + k 2 a 2 v P L {k) from Eq.(53J. In 
other words, we complete the expansion (|53p by substitut- 
ing the terms obtained for the Zeldovich power spectrum 



for all higher orders, which allows us to perform an ex- 
plicit resummation as in the second line of Eq. (|56[) . This is 
clearly a systematic procedure, that can be pushed up to 
any order, and agrees with standard perturbation theory 
up to the order of the substitution. Since the resummed 
Zeldovich power spectrum no longer shows the spurious 
Gaussi an decay of Eg .(1531 but a more physical p ower-law 
decay (|Tavlor fc Hamilton! 119961: IValaeeasI [2007bl l2010al ). 
one may hope that this could improve the agreement with 
results from numerical simulations. 

We show in Fig. [T5] the results we obtain using either 
Ea.([53]l or Ea.([56]>. at redshifts z = 0.35, 1, and 3. As 
pointed out above, we can see that the Gaussian damp- 
ing term of Eq. (1531 leads to a fast decay for the nonlinear 
density power spectrum, which cannot be compensated by 
the 1-halo term. The "modified Zeldovich" expansion (j56|) . 
which only implies a slower power-law decay at high fc, fares 
better as the nonlinear power spectrum keeps growing on 
these scales. However, this is not sufficient to obtain a sat- 
isfactory match with the N-body results. Thus, the com- 
parison with Fig. [10] shows that these Lagrangian-based 
expansions are not competitive with Eulerian schemes, and 
we do not investigate further these approaches in this pa- 
per. We note however that this is rather unfortunate, since 
in principles Lagrangian perturbation theories would seem 
better suited to compute the 2-halo term J2HJ), and more 
importantly they should be bette r suited to computations 
in redshift space (jMatsu bara 2008). Therefore, it remains of 
interest to investigate whether other resummation schemes 
could be developed within the Lagrangian framework, but 
Fig. HU suggests that some important new ingredients are 
required to make significant progress. 



7. Real-space two-point correlation 

Finally, we consider in this section the real-space two-point 
correlation function, £(|x2 — xjj) = (5(xi)#(x2)), that is 
predicted by our model. It is related to the Fourier-space 
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Fig. 20. The real-space two-point correlation function £(x), at redshifts z — 0.35, 1, and 3, as in Fig. [T91 but on larger 
scales. Since the 1-halo contribution is negligible on these scales we only show the linear correlation (black dashed line) 
and the full nonlinear correlation (green solid line). 



power spectrum through 

Ox) = 4tt [°° dkk 2 P(k)^^- (57) 

JO "> x 

= r^ A2(fc) ^£) ) (58) 

so that we obtain £ (x) by integration over k of the nonlinear 
power spectrum (|22p. In particular, from the decomposition 
(|22[) we can write the nonlinear correlation as the sum of 
1-halo and 2-halo contributions, 

=6h(z)+6h(z), (59) 

which are the Fourier transforms of the power spectra (|2"3"|) 
and 

We show our results at redshifts z = 0.35, 1, and 3, 
in Figs. Q1J] and Here we consider our fiducial model 
presented in Sects. [3] and [5j with the 2-halo perturbative 
term P pe rt(fc) given by the direct steepest-descent scheme 
detailed in App. [A] and the 1-halo contribution obtained 
for halos defined by a density contrast <5 = 200 and the 
mass-concentration relation P2")) . 

In agreement with our results for the power spectrum 
shown in Fig. [§1 we can see in Fig. 1191 that our model pro- 
vides a good match to numerical simulations on both small 
and large scales. As is well known, on large scales the two- 
point correlation is governed by the 2-halo contribution £ 2 h, 
whereas on small scales it is governed by the 1-halo contri- 
bution. As for the power spectrum, the transition scales 
where the 2-halo and 1-halo contributions are of the same 
order are more difficult to reproduce, and we again under- 
estimate the power (here the real-space correlation) in this 
range. This is a direct consequence of the same feature ob- 
served in Fourier space in Fig. [HI In particular, we also find 
that the range and amplitude of this discrepancy is larger 
at z = 3 than at z — 0.35. Again, improving theoretical pre- 
dictions on this domain would require pushing perturbation 
schemes to higher orders and/or improving the description 
of halo outer shells. 

We focus in Fig.(2D]on the large scales around the baryon 
acoustic peak. We can check that our model agrees very 
well with the N-body results and reproduces the well-known 



damping of the BAO peak by the nonlinearities of the dy- 
namics (mode-coupling effect). This shows that for the pur- 
pose of BAO studies, the accuracy of our theoretical model 
is certainly sufficient. In particular, on these large scales, 
which are far from the transition region seen in Fig. 1191 the 

1- halo contribution is completely negligible and the two- 
point correlation is set by the perturbative term P per t of the 

2- halo contribution. Thus, it is given by a systematic per- 
turbation theory, and we can check that the simple 1-loop 
direct steepest-descent scheme detailed in App. [A] already 
provides quite satisfactory results. For practical purposes, 
the main improvement that remains to be achieved is to 
extend these results to redshift space, but this is beyond 
the scope of this paper and we leave this point for future 
studies. 

We compare in Fig. [5TJ the correlation functions ob- 
tained with halos defined by a density threshold <5 = 200 
or S = 50. This is the real-space counterpart of Fig. QTJ 
In agreement with Fig. [Til we can see that extending the 
truncation radius of the halo profiles increases somewhat 
the power in the transition range between the 2-halo and 1- 
halo contributions. However, this is not sufficient to obtain 
a perfect match to the numerical simulations, especially at 
z = 3. Again, this suggests that part of the discrepancy is 
due to the incomplete resummation of higher order terms 
of the perturbative contribution P pcrt (fc). 

We do not plot the counterpart of Fig. [20l because on 
these large scales the 1-halo contribution is very small and 
we have checked that the curves obtained with halos defined 
by S = 200 or S — 50 cannot be distinguished in the figure: 
the two curves coincide within the thickness of the curves 
of Fig. [5D1 This also shows that the predictions on these 
large scales are quite robust, since they are independent 
of the underlying halo model and are given by systematic 
perturbation theory. The effect of the 1-halo counterterm 
is also very small so that we do not plot a specific figure to 
compare the curves obtained with and without this coun- 
terterm. Note that the difference between various curves in 
Figs. [12] and [13] was amplified by the fact that we divided 
by a "no- wiggle" linear power spectrum PL S (k). 




Fig. 22. The accuracy of our model and our simulations at redshifts z — 0.35, 1, and 3, for the density power spectrum. 
The green solid line is the relative difference between the model and the simulations, from Eg . ([SO)) . The dashed line that 
grows at low k is the relative statistical error (1451) of the simulations, while the dashed line that grows at high k is the 
relative shot-noise error. 
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Fig. 21. The two-point correlation £(x), as in Fig. \T§\ but 
for halos that are defined either by a density contrast 5 = 
200 (solid curves, identical to those of Fig. lTOl) . or by 6 = 50 
(dot-dashed curves). 

8. Typical accuracy of combined models 

In this last section before conclusion, as a summary of what 
can be achieved from simple analytical models that com- 
bine perturbation theories with halo models, we estimate 
the accuracy reached by our fiducial model. Thus, we plot 
in Fig. [52] the relative difference between our model, de- 
scribed in Sects. [3] and [5J and the results from the N-body 
simulations presented in Sect. |4j 

AP = |Pmodol(fc)-ifr-body(fc)| (6Q) 
P -PN-body(fc) 

We also show the statistical error due to the finite num- 
ber of modes in the simulation box, given by Eq. (l45p . and 
the shot-noise error given by AP^ot- noise = L^ oyi /N . Since 
the number of modes within a bin of fixed size Ak scales 



as k 2 Ak, the relative statistical error decreases at higher 
k as 1/k, see Eq. (|4"5")) . The sudden jumps are due to the 
folding procedure, see Sect. 14.21 and Fig. [7] On the other 
hand, the relative shot noise grows as 1/P(k) and domi- 
nates at high k. It seems that the simple approximation 
AP s hot-noise = L^ ox /N overestimates somewhat the inac- 
curacy of the simulations, however we do not look for a bet- 
ter estimate here, as this is already sufficient to understand 
the high-fc part of Fig. [2D and to mark the wavenumber 
where shot noise becomes dominant. 

We can see that on the large scales associated with the 
BAO oscillations, 0.03 < k < 0.3/iMpc -1 , the analytical 
approach is competitive with the numerical simulations. 
It typically gives an accuracy of 1%, which however wors- 
ens somewhat at lower redshift to reach 2% at z = 0.35. 
However, a part of AP/P shown in Fig. |2"21 arises from the 
inaccuracies of the simulations themselves, so that the true 
value of relative error between the analytical model and the 
exact nonlinear power spectrum may be slightly smaller. Of 
course, on even larger scales, k < 0.03/iMpc -1 , the analyt- 
ical model becomes exact as it agrees with linear theory 
and the relative difference AP/P is solely due to the in- 
accuracies of the simulation results (as can be checked by 
comparison with the relative statistical error of the simula- 
tions). In the 1-halo region, 3 < k < 100/iMpc^ 1 , we reach 
an accuracy on the order of 10%, and even better (4%) in 
the region where 10 < A 2 (fc) < 100. Moreover, as we have 
seen in Sect. [5J the predictions of our model are rather ro- 
bust in this range, since the dependence on the details of 
the underlying halo model is quite weak. At very high k the 
numerical simulations are strongly affected by shot noise, 
so that the curve AP/ P mostly measures this source of er- 
ror and the accuracy of the model may be better than what 
Fig. suggests. 

In agreement with the discussion in Sect. [5l a salient 
feature in Fig. l22l is the peak in the transition range between 
the 2-halo and 1-halo contributions, at k ~ 1/iMpc -1 . This 
is due to the underestimation noticed in Fig. |9l which can 
reach 30% around the peak. This feature shifts to higher 
k and somewhat broadens at higher redshift, in agreement 
with the discussion in Sect. \b[ 

Next, we show in Fig. [531 the relative accuracy of our 
model and our simulations for the real-space two-point cor- 
relation function. As in the Fourier-space figure we 




Fig. 23. The accuracy of our model and our simulations at redshifts z — 0.35, 1, and 3, for the density two-point 
correlation. The green solid line is the relative difference between the model and the simulations, from Ea. ljim) . The 
dashed line that grows at large x is the relative statistical error (|46|) of the simulations, while the dashed line that grows 
at small x is the relative shot-noise error. 



plot the relative difference between our model, presented in 
Sect. [7] (for halos defined by 8 = 200), and our simulations, 

A£ . . _ |£model(gO - ^N-bodyQc)! 
/■ \ X ) ~ (■ ( \ ■ x *-) 

5 4N-bod y (a;j 

We also show the statistical error, given by Ea. (l46l) . 
and the shot-noise error given by A£ sno t- noise = 
[L^ ox /(47ra; 3 /3)]/iV. As for the power spectrum, the rela- 
tive statistical error grows on large scales, as oc x, because 
of the smaller number of modes, and the sudden jumps are 
due to the folding procedure, see Sect. 14.21 and Fig. [JJ On 
the other hand, the relative shot noise grows as l/(x 3 £,(x)) 
on small scales and dominates at low x. 

We can check that Fig. [23] is consistent with Fig. [2U 
On large scales, x > 10/i _1 Mpc, the analytical approach 
is again competitive with the numerical simulations and it 
reaches an accuracy of 1%. On very large scales the analyti- 
cal model converges to the linear two-point correlation and 
becomes exact, and the relative difference A£/£ is solely 
due to the inaccuracies of the simulations. On small scales 
dominated by the 1-halo contribution, x < 0.1/i _1 Mpc, we 
reach an accuracy of 10%, and even better at z > 1. On 
very small scales the difference A£/£ mostly measures the 
level of shot noise in the simulations, so that the accuracy 
of the model could be better than shown in the figure. 

Again, the most difficult region to reproduce is the 
transition between the 2-halo and 1-halo contributions, 
x ~ l/i _1 Mpc. The underestimation seen in Fig. leads 
to a peak for the relative error of the analytical model that 
can reach about 30% and is somewhat greater at higher 
rcdshift. Of course, the boundaries of these various domains 
shift to smaller scales at higher redshifts. 

Thus, it appears that analytical models such as the one 
studied in this paper can already provide reasonably good 
estimates for the density power spectrum and two-point 
correlation function over a large range of scales. We have 
not scanned all possible models in this paper (although we 
have investigated several variants in Sect.|5]), but we expect 
similar approaches, based on a combination of perturbation 
theory and halo model, to yield similar levels of accuracy, 
for current resummation schemes and halo models. At low 
and high k, and at small and large x, the level of accuracy 
is already rather satisfactory, and may be improved in a 



straightforward and systematic manner by pushing pertur- 
bation theories to higher orders and measuring halo profiles 
down to smaller scales in simulations (and possibly taking 
into account additional effects such as substructures and 
baryon impact on gravitational collapse). An advantage of 
such models, especially on large scales, it to provide robust 
predictions that can be applied to a variety of initial con- 
ditions or cosmological parameters. This is obvious on the 
large scales described by systematic perturbation theories, 
but this remains true to some extent at high k and small 
x as the power spectrum is expressed in terms of more el- 
ementary quantities that satisfy physical constraints (such 
as the normalization or scaling of the mass function, and 
reasonable halo profiles) . This ensures at the very least that 
the power spectra and two-point correlations obtained in 
this fashion make sense and scale in the appropriate fash- 
ion, even for (reasonable) cases where they have not been 
tested. 



What may prove to be the most difficult regime for sys- 
tematic improvement is the transition range. Although we 
have not investigated such an approach here, it is certainly 
possible to build fitting formulas that improve the agree- 
ment with simulations on this range. Thus, one may choose 
some interpolation formula between the low-fc and high- 
k domains and obtain a good match by tuning some free 
parameters that govern the curvature of the interpolation 
curve in this regime. However, one must take some care not 
to spoil the good agreement on BAO scales and to keep a 
percent accuracy there. On the other hand, if we insist on 
avoiding such a direct interpolation procedure, so as to de- 
rive the power spectrum from more elementary quantities 
(such as perturbative expansions and halo profiles) with- 
out further modifications, it is not clear how much progress 
can be made on this range. Indeed, the description of the 
density field in terms of halos is by itself an approximate 
picture, especially for the moderate-density regions of space 
associated with outer shells or filaments that have already 
experienced shell crossing but have not fully virialized and 
are not enclosed within well-identified objects. This may 
prevent a perfect matching on the transition scales, unless 
it is somehow enforced a posteriori. 
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9. Conclusion 

In this article we have explained how to build unified models 
to predict the matter density power spectrum, and the real- 
space two-point correlation function, by combining pertur- 
bation theories with halo models. First, we have shown 
how a Lagrangian point of view allows us to re-interpret 
the decomposition of the power spectrum into 2-halo and 
1-halo contributions. This provides a convenient route to 
this splitting, which automatically ensures the conserva- 
tion of mass and gives an explicit relationship with the dy- 
namics of the system, through the direct Lagrangian map 
q M> x. We have also emphasized the relationship between 
this decomposition and the separation into perturbative 
and non-perturbative terms. This explicitly shows how one 
can achieve consistency with perturbation theory through 
the 2-halo term, as the 1-halo term is fully non-perturbative 
and does not contribute to any order of perturbation the- 
ory. 

Next, we have shown that the 1-halo contribution con- 
tains a counterterm, which was missed in previous studies, 
that gives the expected low- A; tail Pm(k) oc k 2 associated 
with the conservation of matter. This ensures that the 1- 
halo term becomes subdominant on large scales, as linear 
CDM power spectra typically behave as Ph{k) oc k at low 
k. On the other hand, we have explained why the fc 4 -tail 
associated with momentum conservation is not recovered 
because of the approximations that we use (such as instan- 
taneous virialization and loss of memory), but this is not a 
problem for practical purposes. If needed, in future works 
aiming at greater accuracy that include higher orders of 
perturbation theory, it may be possible to enforce the k - 
tail by modifying this counterterm. Then, we have pointed 
out that the 2-halo contribution, and more precisely its 
term P pe rt(k) that should be consistent with perturbation 
theory, cannot be computed by truncating standard pertur- 
bation theory at a finite order. Indeed, higher order terms of 
this perturbative expansion grow increasingly fast at high k 
(with large cancellations between various orders), so that an 
abrupt truncation at a given order leads to a 2-halo contri- 
bution that is non-negligible at high k (and would even be 
dominant on small scales if we go to two loops or higher or- 
ders), which is not physical and prevents a good agreement 
with numerical simulations. The remedy is to use resumma- 
tion schemes that agree with standard perturbation theory, 
up to the required order, while remaining well-behaved at 
high k (typically close to or smaller than the linear power 
spectrum). In addition to the superior accuracy of such 
schemes on large scales, this is a second important moti- 
vation to develop such perturbative approaches. We have 
described in particular a simple implementation, based on 
the direct steepest-descent resummation at 1-loop order, 
which satisfies this property while being fast to compute, 
thanks to its factorized form. In this fashion, both the 1- 
halo and 2-halo contributions remain well-behaved beyond 
the domain where they dominate. This allows us to build 
a meaningful model from their combination, which can de- 
scribe all scales and regimes. 

We have compared a simple implementation of this 
model to N-body simulations. We can reach an accuracy 
of 1% on quasi-linear scales, associated with baryon acous- 
tic oscillations, and 10% on highly nonlinear scales domi- 
nated by the 1-halo contribution. The agreement on large 
scale is actually better than the one obtained with stan- 



dard 1-loop perturbation theory. Since it is based on sys- 
tematic perturbation theory and contains no free parameter 
(except through the prefactor i<2H, which however is very 
close to unity on these scales), this gives a robust and re- 
liable prediction that also improves over simple fitting for- 
mulas. On the other hand, the good agreement on highly 
nonlinear scales, A 2 (fc) > 200, depends on the properties 
of the halos included in the underlying halo model, and 
more specifically on the mass-concentration relation c(M). 
We have given a simple model that provides a reasonable 
match to the N-body results, but this regime is clearly less 
reliable than the quasi-linear regime. In particular, to ap- 
ply this framework to very high redshifts one should use 
a better constrained prescription c(M) or a physically mo- 
tivated model. Nevertheless, in this high-fc regime where 
200 < A 2 (fc) < 1000, we find that using various formulas 
for c(M) that have been proposed in the literature already 
improves over simple fitting formulas for P(k), in spite of 
the scatter of their predictions. Indeed, on these scales the 
dependence on the details of the model (the low-mass tail 
of c(M)) is not too large yet. Moreover, on somewhat less 
nonlinear scales where 10 < A 2 < 200, the predictions ap- 
pear even more robust and also provide a good match to 
simulations. 

The transition range where the 2-halo and 1-halo contri- 
butions are of the same order, roughly where 1 < A 2 < 10, 
proves to be the mos t difficult to reprod uce. In agreement 
with previous works (jGiocoli et al.l [20ltf). we find that the 
model underestimates the power spectrum on these scales, 
and the inaccuracy is somewhat larger at higher z. Some 
of the discrepancy may arise from the truncation of halo 
profiles at a density contrast 6 — 200, and we have shown 
that going to larger radii, defined by S = 50, gives a slightly 
better match to the simulations, but this is not sufficient. 
Another route to improve the model would be to truncate 
the resummation scheme (involved in the 2-halo contribu- 
tion) at a higher order, that is to include exactly all pertur- 
bative terms up to a higher loop order, while higher orders 
are only partial ly resummed. As seen in previous works 
(|Valageasll2010al) , the potential for improvement is greatest 
at higher z for CDM power spectra, since the range where 
perturbative expansions are relevant is somewhat broader 
and one can push to higher orders before non-perturbative 
terms become dominant. This would explain why the un- 
derestimate of the power spectrum appears to be more sig- 
nificant at z = 3 than at z — 0.35. 

However, it is likely that even a complete resummation 
of the perturbative term would not provide a perfect match, 
because of the approximations involved by the splitting into 
2-halo and 1-halo contributions. Indeed, even though the 
decomposition (|2"2"j). with Eas. (|2"51) - (|24D . can be made ex- 
act if we have a well-defined criterion to define halos, the 
computation of the 2-halo term involves some approxima- 
tion as we replace the conditional average (-.)2H by the 
perturbative average (..} pe rt- Going beyond this approxi- 
mation requires going beyond perturbation theory (which 
does not know about halos). On the other hand, the 1-halo 
contribution also involves approximations that are difficult 
to improve, such as the instantaneous virialization, that is, 
the fact that particles are located at random in a halo, in- 
dependently of their initial location. Nevertheless, this is a 
range of scales where there is still room for improvement. 

Next, we have investigated the dependence of such pre- 
dictions on the details of the model. We have shown that 
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it is important to take into account the 1-halo counterterm 
in order to get a good accuracy on large scales, k ~ 0.1 
to 0.3/iMpc -1 . Otherwise, the 1-halo contribution does not 
decrease fast enough at low k and spoils the good match 
with simulations (and would even become dominant on very 
large scales). On the perturbative side, we have shown that 
standard perturbation theory cannot be used to build uni- 
fied model, because it yields higher order terms that grow 
increasingly fast at high k and prevent a good match to nu- 
merical simulations in the nonlinear regime. This cannot be 
compensated by details of the underlying model and would 
require at best the addition of some ad-hoc damping prcf- 
actor. Another perturbative scheme than the 1-loop direct 
steepest-descent resummation we have focused on in this 
paper, based on a more realistic "decaying response" func- 
tion, gives similar results to the former, as it also provides 
a well-behaved 2-halo contribution at high k. It even im- 
proves somewhat the match with simulations on transition 
scales, but this is likely to be a coincidence as this comes at 
the price of an overestimate on the larger scales probed by 
BAO. Next, we have shown that current Lagrangian resum- 
mation schemes, which a priori could have been expected 
to be best suited to our purposes, actually fare much worse 
than Eulerian resummation schemes. Indeed, although they 
yield a well-behaved perturbative term at high k, the asso- 
ciated damping is so large that it leads to a significant un- 
derestimate for the power spectrum already on BAO scales, 
k ~ 0.1 to 0.3/iMpc -1 . However, if it is possible to devise 
new Lagrangian schemes that do not suffer from this ef- 
fect they may provide a promising route to improve unified 
models for the power spectrum. 

Finally, we have shown that our model also provides 
reliable predictions for the density two-point correlation 
function. As for the power spectrum, it describes rather 
well both small and large scales, but underestimates some- 
what the correlation in the transition range between the 2- 
halo and 1-halo contributions. On very large scales, around 
the BAO peak, we obtain a very good agreement with nu- 
merical simulations, so that the simple 1-loop resummation 
scheme used in this work is probably sufficient to describe 
this feature. 

Thus, we have shown that it is possible to build suc- 
cessful unified models for the density power spectrum and 
two-point correlation, that combine systematic and reliable 
perturbation theories on large scales with phenomenological 
halo models on small scales. Two issues where there is still 
room for improvement are the underestimate of the power 
spectrum and correlation function on transition scales and 
the building of accurate Lagrangian resummation schemes 
(but this is a more theoretical goal since for practical pur- 
poses Eulerian schemes are sufficient). A third issue is the 
modeling of the underlying virialized halos, which plays a 
key role at very high k. There, many points can be im- 
proved, such as taking into account substructures, devia- 
tions from spherical profiles, or the effect of baryons on 
the halo density profile. In particular, one should be cau- 
tious when using such models beyond the range where they 
have been tested (A 2 (fc) < 1000 and z < 3). There, the 
main uncertainty comes from the low-mass tail of the mass- 
concentration relation c(M), and it is probably safer to use 
physically motivated models for c(M) (even though they 
may not be very accurate), or at least to check that the 
formula used for c(M) is still reasonable. 



The model presented in this paper could be extended 
to higher order statistics, such as the bispectrum or three- 
point correlation function. It should also be possible to han- 
dle the case of non-Gaussian initial conditions, at least as 
long as the halo density profiles are not strongly modified 
(or this would require a separate study to evaluate this 
effect). Finally, to allow a direct comparison with galaxy 
surveys, it would be desirable to extend the model to red- 
shift space. This is likely to be a more involved task, that 
we leave for future works. 
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Appendix A: Computation of P pe rt{k) with the 
"direct steepest-descent method" 

We briefly recall here how the perturbative term P pe rt(k) of 
the 2-halo contribution (f2"7| is com puted within the " direct 
steepest-descent method", see also IValaeeasI (|2007aL [2008) 
for details. First, to simultaneously describe the density and 
velocity fields it is convenient to define the two- component 
vector V>(k, rj) as (jCrocce fc Scoccimarrdl20"06b[ ). 



6(k, rj) 
-0(k, v )/(af) 



(A.l) 



where 9 is the divergence of the velocity field, 9 = V • v, 
i] = \nD(z) where D(z) is the linear growth factor (nor- 
malized by D(0) = 1 today), d = da/dt where a(t) is the 
scale factor, and / = dlnZ?/dlna. Then, introducing the 
coordinate x — (k, 77, i), where i — 1 or 2 is the discrete in- 
dex of th e two-componen t vectors, the equations of motion 
write as dWagcas 20072) 



0(x,x') ■ ip(x') = K s (x;xi,x 2 ) ■ ip(xi)ip(x 2 ), 
where the matrix O reads as 



(A.2) 



<D{x,x') 



8 D {k-k!)5 D {n-n') 
( 

dn 
30„ 



X 



2/ 2 



dn 



-1 



1 



(A.3) 



whereas the symmetric vertex K s (x] xi,x 2 ) = K s (x; x%,xx) 
writes as 



K s (x;xi,x 2 ) 
with 

7i s ; i :2 ( k i> k 2) 

7l;2,l( k l> k 2) 
7 2;2 , 2 ( k l, k 2) 



5 D (ki + k 2 - k) 5 D (r)i 

X 7k, l2 ( k l' k 2), 

(ki + k 2 ) ■ k 2 

2fc| 

(ki + k 2 ) ■ ki 

2fc2 

|k!+k 2 | 2 (ki -k 2 ) 

2 }\>-y k<2 



rf) S D (ri2 - rf) 

(A.4) 



(A.5) 
(A.6) 
(A.7) 



P. Valageas and T. Nishimichi: Combining perturbation theories with halo models 



27 



and zero otherwise (|Crocce &: Scoccimarrol l2006bf) . In ad- 
dition to the usual two-point correlation function, 



C(xi,x 2 ) = (ip(xi)ip(x 2 )), 



(A.8) 



resummation schemes often involve the response function 
(or propagator) , defined as the functional derivative of the 
nonlinear field ip with respect to an infinitesimal "noise" £ 
that would be added to the right hand side of the equation 
of motion (|A.2|) . 



Mxi) 



(A.9) 



As described in IValageasI (|2007al I2008D . the nonlin- 
ear correlation and response functions obey the exact 
Schwinger-Dyson equations 

0(x,z)-C(z,y) = Y.(x,z) ■ C(z,y) + U(x,z) ■ R T (z,y), 

(A.IO) 

0(x, z) ■ R(z, y) = 5 D (x-y) + Efo z) ■ R(z, y), (A.ll) 

where E(x,j/) and H(x,y) are "self-energy" terms, that can 
be written as diagrammatic expansions over the linear two- 
point functions Cl and Rl- Equation (lA.lOp can be inte- 
grated as 



C = Rx C L (r/i) x R T + R ■ II • R 1 



(A.12) 



where Cl{i1i) is the linear correlation function at an early 
time 7/1, with 7/1 — oo, and the first product "xCix" 
does not contain any integration over time. 

Then, at 1-loop order the direct steepest-descent scheme 
amounts to use for E and II their lowest order terms, over 
powers of Cl and Rl- This means that Eq. (|A.12[) for the 
correlation C reads as the first line in Fig. [TJ while a sim- 
ilar diag rammatic expre ssion holds for the response R (see 
Fig.8 in IValageasI (|2008l) h In particular, solving Eg. ([AT f ft 
as a pcrturbative series over powers of E gives a diagram- 
matic expansion for the response R in terms of "bubble" 
diagrams, and substituting into Eq. (|A.12[) this gives the 
series of diagrams shown in the second equality of Fig. [T] 
In practice, one does not compute the two-point function 
C(x, y), whence the density power spectrum, from these di- 
agrammatic series, but from the integro-differential equa- 
tion (|A.I1[) and the explicit expression (|A.12[) . 

First, using in the following the approximation 
^m// 2 — 1, the linear correlation and response functions 
read as 



C L (si , x 2 ) = e" 1 +" 2 P L0 (h ) S D (ki + k 2 



1 I 
1 I 



(A.13) 



where PLo(k) is the linear density power spectrum today 
(when 77 = 0), and 

Rl(xi,x 2 ) = 0(771 - r] 2 ) 5 D (ki - k 2 ) 

e ,11 - TI2 R+ + e~ 3 ^-^/ 2 Ro] , (A.I4) 



with 



+ _ / 3/5 2/5 \ _ / 2/5 -2/5 \ 

U ° ~ \ 3/5 2/5 J and U ° ~ { -3/5 3/5 J ' 



15) 



In Eq. (|A.14p the Heaviside factor 0(771 — 772) expresses the 
causality of the response funct ion. Then, the " self-energy" 
terms read at lowest order as (jVala acas 2007a|) 

T, (xi,x 2 ) = 0(771 - 772) S D (ki - k 2 ) 

e^E+^i) + e -m/2+5 % /2 s - (fci) l i (A lg) 



with 



J 0;iii 2 



(k) 



4 J2 J d q7l i;Jlj2 (k-q,q) 



X 



7f i;i2 ,/ 2 (k,-q)P L0 ( g )i?± iiii , 



and 

Ho (3:1,0:2) 
with 



fe(ki+k 2 )e 2 " 1 + 2 " 2 n (fei) 



(A.17) 
(A.18) 



n, 



0;iii 2 



(k) = 2 ]T J dq7f i;Jlj2 (q,k-q) 

jitehh 

xil;i uh ("Ob -k + q) PM P L o(\k - q|). (A.19) 
These "self-energy" terms also satisfy the symmetries 

En 



So (*) = 



J 0;22 
-E 



-E+ 

^0:12 



E+ 

0:21 ^0:11 



U . lJ {k)=U . Jl (k). (A.20) 



Thanks to the simple time dependence of Eo(xi, x 2 ), which 
arises from the fact that we have expanded E over the linear 
two-point functions Cl and Rl and only kept the lowest 
order term, it is possible to eliminate the integral terms in 
the Schwinger-Dyson equation (|A.11[) by taking successive 
derivatives with respect to 771 . This yields 



<9 3 i?i 3 d 2 Ri 



l 2 dRi 3 dR 2 



P 2m 



2 dr/ 2 dril drji 2 $771 



Ri = 



3?7i 



^0;22 



9^11^1 + 9^0^12-^2 



d 3 R 2 3d 2 R! d 2 R 2 9 5i?i 



7dR 2 



dm 

c 2r,i 



2 <9t7 2 

(Km ^ 



J 0;22 



d-q 2 4 <9?7i 

, dR 2 
dm ' 2 



4 9771 

-Ept 21 iii + -E(}: 22 i? 2 



(A.21) 



\Ri ~ ^2 - 



(A.22) 



where the pair (i?i, i? 2 ) corresponds either to (i?n, i? 2 i) or 
(i?i2, i? 22 )- These equations apply to 771 > 772 and the initial 
conditions at 771 = 772 are 



X 1 



Ri R 2 
1 



dRi 8R 2 d 2 B 

dr/i drji dr/'f drj^ 



2 R2 



3 3 



+ e ^(E+ 11 + E+ 22 ) -|) 



for the pair (i?n,i? 2 i), and 



r = 011 



1 1 

2 4 



o2r,i 



(^0:11 



J 0;22/ 



(A.23) 
(A.24) 

(A.25) 



for the pair (Ri 2 ,R 22 ). Thus, we obtain the response func- 
tion by solving the differential equations (|A.21[) - (|A.22|) . 
with the initial conditions (IA.24|l - (jA.25|l . using an ordinary 
Runge-Kutta algorithm. A great simplification is that there 
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are no longer integral terms over past times in the right 
hand side of Eos. (|A.2ip - (|A. 22ft . so that numerical compu- 
tations are very fast. 

Finally, it is possible to take advantage of the factorized 
time dependence of Ilo in Eq. (|A.18|) to factorize Eq. (|A.12l) . 
To compute the first product in Eg. (jA. 12ft we define the 
two-component vector 

2 

V\{k, V) = J2 R ^ 1> »») em VPMk), (A.26) 
i=i 

while for the second product we diagonalize the symmetric 
matrix Ilo as 

II (fc) =p+ 7T+ -Tr +T + P - it- -7r- T , (A.27) 

where p^ are the eigenvalues of ITo and tt^ its normalized 
eigenvectors, and we define the two-component vectors, 

rv 

(p ± (k,rj)= dr]'e 2ri R{k; 77,77') ■w ± (k). (A.28) 
J — 00 

Then, Eq. (|A.12l) writes as 
C(fc;77i,7j 2 ) = </(rn)-<P I (V2) T 

+ P +tp+( m ) ■tp+im f + p -^-( m ) -tp-im f- (A.29) 

In particular, the nonlinear power spectrum is given by the 
equal-time upper-left component of the correlation matrix 

C, 

P(k,r)) = Cn(k;v,v)- (A.30) 

The factorized form (|A.29[) . which avoids to compute two- 
dimensional integrals over past times as in generic expres- 
sions of Ea. (|A.12|) . thanks to the factorized time depen- 
dence of n in Eq. (|A.18p . again allows fast numerical com- 
putations. 

References 

Abramowitz, M. & Stegun, I. A. 1970, Handbook of Mathematical 

Functions (New York: Dover) 
Bernardeau, F., Colombi, S., Gaztanaga, E., & Scoccimarro, R. 2002, 

Phys. Rep., 367, 1 
Bernardeau, F. & Valageas, P. 2008, Phys. Rev. D, 78, 083503 
Bernardeau, F. & Valageas, P. 2010a, Phys. Rev. D, 81, 043516 
Bernardeau, F. & Valageas, P. 2010b, Phys. Rev. E 
Bouchet, F. R., Colombi, S., Hivon, E., & Juszkiewicz, R. 1995, 

Astron. Astrophys., 296, 575 
Buchert, T. 1994, Mon. Not. R. Astron. Soc, 267, 811 
Carlson, J., White, M., & Padmanabhan, N. 2009, Phys. Rev. D, 80, 

043531 

Colombi, S., Jaffe, A., Novikov, D., & Pichon, C. 2009, Mon. Not. R. 

Astron. Soc, 393, 511 
Cooray, A. & Sheth, R. 2002, Phys. Rep., 372, 1 

Crocce, M., S. Pueblas, S., & Scoccimarro, R. 2006, Mon. Not. R. 

Astron. Soc, 373, 369 
Crocce, M. & Scoccimarro, R. 2006a, Phys. Rev. D, 73, 063520 
Crocce, M. & Scoccimarro, R. 2006b, Phys. Rev. D, 73, 063519 
Crocce, M. & Scoccimarro, R. 2008, Phys. Rev. D, 77, 023533 
Desjacques, V. & Seljak, U. 2010, Classical and Quantum Gravity, 27, 

124011 

Dolag, K., Bartelmann, M., & et al., F. P. 2004, Astron. Astrophys., 
416, 853 

Duffy, A. R., Schaye, J., Kay, S. T., & Vecchia, C. D. 2008, Mon. Not. 

R. Astron. Soc, 390, L64 
Eisenstein, D. J. & Hu, W. 1999, Astrophys. J., 511, 5 
Eisenstein, D. J., Hu, W., & Tegmark, M. 1998, Astrophys. J. Lett., 

504, 57 



Eisenstein, D. J., Zehavi, I., & et al., D. W. H. 2005, Astrophys. J., 
633, 560 

Feldman, H. A., Kaiser, N., & Peacock, J. A. 1994, Astrophys. J., 426, 
23 

Giocoli, C, Bartelmann, M., Sheth, R. K., & Cacciato, M. 2010, Mon. 

Not. R. Astron. Soc, 
Goroff, M. H., Grinstein, B., Rey, S.-J., & Wise, M. B. 1986, 

Astrophys. J., 311, 6 
Gurbatov, S. N., Saichev, A. I., & Shandarin, S. F. 1989, Mon. Not. 

R. Astron. Soc, 236, 385 
Hamilton, A. J. S. 1992, Astrophys. J. Lett., 385, 5 
Hamilton, A. J. S., Kumar, P., Lu, E., & Matthews, A. 1991, 

Astrophys. J. Lett., 374, 1 
Hockney, R. W. & Eastwood, J. W. 1981, Computer Simulation Using 

Particles (New York: McGraw-Hill) 
Huffenberger, K. M. & Seljak, U. 2003, Mon. Not. R. Astron. Soc, 

340, 1199 

Komatsu, E., Dunkley, J., & et al., M. R. N. 2009, Astrophys. J. 
Suppl., 180, 330 

Lewis, A., Challinor, A., & Lasenby, A. 2000, Astrophys. J., 538, 473 
Liguori, M., Sefusatti, E., Fergusson, J. R., & Shellard, E. P. S. 2010, 
arXiv:1001.4707 

Massey, R., Rhodes, J., & et al., A. L. 2007, Astrophys. J. Supp., 172, 
239 

Matarrese, S. & Pietroni, M. 2007, JCAP, 6, 26 
Matsubara, T. 2008, Phys. Rev. D, 77, 063530 

Mo, H. J. & White, S. D. M. 1996, Mon. Not. R. Astron. Soc, 282, 
347 

Munshi, D., Valageas, P., van Waerbeke, L., & Heavens, A. 2008, 

Phys. Rep., 462, 67 
Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, Astrophys. J., 

490, 493 

Peacock, J. A. & Dodds, S. J. 1996, Mon. Not. R. Astron. Soc, 280, 
L19 

Peebles, P. J. E. 1974, Astron. Astrophys., 32, 391 

Peebles, P. J. E. 1980, The large scale structure of the universe 

(Princeton: Princeton university press) 
Peebles, P. J. E. 1982, Astrophys. J. Lett., 263, 1 
Percival, W. J. & White, M. 2009, Mon. Not. R. Astron, 393, 297 
Pietroni, M. 2008, JCAP, 10, 36 

Press, W. & Schechter, P. 1974, Astrophys. J., 187, 425 
Saito, S., Takada, M., & Taruya, A. 2010, arXiv:1006.4845 
Scherrer, R. J. & Bertschinger, E. 1991, Astrophys. J., 381, 349 
Schneider, P. & Bartelmann, M. 1995, Mon. Not. R. Astron. Soc, 
273, 475 

Scoccimarro, R. 1998, Mon. Not. R. Astron. Soc, 299, 1097 
Scoccimarro, R., Sheth, R. K., Hui, L., &; Jain, B. 2001, Astrophys. 
J., 546, 20 

Seljak, U. 2000, Mon. Not. R. Astron. Soc, 318, 203 
Sheth, R. K., Mo, H. J., & Tormen, G. 2001, Mon. Not. R. Astron. 
Soc, 323, 1 

Sheth, R. K. & Tormen, G. 1999, Mon. Not. R. Astron. Soc, 308, 119 
Smith, R. E., Peacock, J. A., & et al., A. J. 2003, Mon. Not. R. Astron. 
Soc, 341, 1311 

Smith, R. E., Scoccimarro, R., & Sheth, R. K. 2007, Phys. Rev. D, 
75, 063512 

Springel, V. 2005, Mon. Not. R. Astron. Soc, 364, 1105 
Swanson, M. E. C, Percival, W. J., & Lahav, O. 2010, arXiv:1006.2825 
Takahashi, R., Yoshida, N., & et al., M. T. 2009, Astrophys. J., 700, 
479 

Taruya, A. & Hiramatsu, T. 2008, Astrophys. J., 674, 617 

Taruya, A., Nishimichi, T., Saito, S., & Hiramatsu, T. 2009, Phys. 

Rev. D, 80, 123503 
Taylor, A. N. & Hamilton, A. J. S. 1996, Mon. Not. R. Astron. Soc, 

282, 767 

Tinker, J. L., Robertson, B. E., & et al., A. V. K. 2010, Astrophys. 
J., 724, 878 

Valageas, P. 2002a, Astron. Astroph., 382, 412 
Valageas, P. 2002b, Astron. Astroph., 382, 450 
Valageas, P. 2004, Astron. Astroph., 421, 23 
Valageas, P. 2007a, Astron. Astroph., 465, 725 
Valageas, P. 2007b, Astron. Astroph., 476, 31 
Valageas, P. 2008, Astron. Astroph., 484, 79 
Valageas, P. 2009a, Astron. Astrophys., 508, 93 
Valageas, P. 2009b, J. Stat. Phys., 137, 729 
Valageas, P. 2009c, J. Stat. Phys., 134, 589 
Valageas, P. 2010a, arXiv: 1009.0106 
Valageas, P. 2010b, arXiv: 1009. 1131 



P. Valageas and T. Nishimichi: Combining perturbation theories with halo models 



29 



Valageas, P. & Bernardeau, F. 2010, arXiv: 1009. 1974 

Vergassola, M., Dubrulle, B., Frisch, U., & Noullcz, A. 1994, Astron. 

Astrophys., 289, 325 
Zeldovich, Y. B. 1970, Astron. Astrophys., 5, 84 



